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1.14.1  INTRODUCTION 

Design  of  composite  aggregates  for  a  specific 
application  involves  several  physical  and 
geometrical  parameters  related  to  the  micro¬ 


structure,  and  depends  on  the  type  of  service 
loads  and  the  application  environment.  More¬ 
over,  fabrication  methods  and  parameters 
cause  significant  local  stresses  and  deformation 
that  may  affect  the  performance  of  the  material. 
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Under  these  circumstances,  it  is  essential  to 
develop  micromechanical  theories  which  evalu¬ 
ate  the  local  fields  and  predict  the  overall 
response  under  combined  thermal  and  mechan¬ 
ical  loads.  The  basic  elements  of  these  theories 
are  geometrical  modeling  of  the  microstruc¬ 
tures  and  local  interactions,  and  constitutive 
modeling  of  the  homogeneous  phases.  An  es¬ 
sential  requirement  of  the  latter  is  to  include 
inelastic  deformation  to  model  composite  sys¬ 
tems  which  exhibit  nonlinear  response  under 
thermal  and  mechanical  service  loads. 

Since  the  early  1960s,  micromechanics  of 
composite  materials  has  attracted  many  re¬ 
searchers.  This  led  to  significant  theoretical 
developments  for  prediction  of  elastic  and  in¬ 
elastic  constitutive  response,  and  motivated  ex¬ 
perimental  validation.  The  purpose  of  this 
chapter  is  to  summarize  these  developments, 
and  illustrate  their  application  in  predicting 
the  overall  response  under  thermomechanical 
loads.  Although  treatment  of  the  subject  in  this 
chapter  covers  two-phase  particulate  and  fi¬ 
brous  materials,  the  focus  in  applications  will 
be  on  fibrous  composites  and  laminates. 

We  begin  by  examining  elastic  composites 
and  provide  general  relations  among  the  local 
and  overall  fields.  This  is  followed  by  descrip¬ 
tion  of  two  classes  of  models,  which  associate 
these  relations  with  an  idealized  microgeometry 
of  the  composite.  In  one  class  of  models,  the 
microgeometry  is  not  known  in  detail,  but  in¬ 
teraction  of  the  phases  is  considered  by  the 
solution  of  certain  inclusion  problems.  In  the 
other  class  of  models,  the  microstructure  is 
idealized  as  a  periodic  dispersion  of  the  reinfor¬ 
cing  phase  into  the  matrix,  and  a  unit  cell  is 
analyzed  under  the  applied  loads.  The  next  two 
sections  are  devoted  to  application  of  these 
models,  when  rate-independent  and  rate-depen¬ 
dent  flow  theories  describe  constitutive  beha¬ 
vior  of  the  phases.  We  close  with  a  section  on 
inelastic  laminates  where  overall  constitutive 
equations,  which  utilize  micromechanical  mod¬ 
els  of  the  individual  plies,  are  developed. 

The  notation  used  here  are  symbolic,  where 
symmetric  second-order  tensors  are  written  as 
(6x1)  matrices  and  denoted  by  boldface  lower 
case  letters,  and  symmetric  fourth-order  tensors 
are  written  as  (6  x  6)  matrices  and  denoted  by 
boldface  upper  case  letters.  Connections  with 
tensor  notation  are  easily  established.  For  ex¬ 
ample,  the  stress  tensor,  a/y,  and  strain  tensor 
s ij\  with  the  symmetry  =  oy-  and  e ~  tjh  are 
written  in  a  matrix  form  as  a  =  [a{l,  ct22,  033, 
°23i  a3b  Gnh  and  £  =  hi,  £221  £33*  2£23,  283j, 
2e  1 2] .  Similarly,  fourth-order  tensors  having  at 
least  the  symmetries  Aijk[  =  Ajikl  —  Aij!k  are  re¬ 
duced  to  (6  x  6)  matrices  A,  such  that 
AA"1  =  A”1  A  =  I,  the  identity  matrix. 


1.14.2  MICROMECHANICS  OF  ELASTIC 
COMPOSITES 

In  this  section,  we  outline  the  general  proce¬ 
dure  for  evaluation  of  the  local  stresses  and 
strains,  as  well  as  the  overall  properties  of  a 
composite  medium  consisting  of  two  distinct, 
elastically  anisotropic  constituents.  The  consti¬ 
tuents  are  assumed  to  be  fully  bonded  at  a 
defined  interface,  and  free  of  voids  and  cracks, 
initially  and  under  subsequent  deformation. 
Loading  consists  of  overall  uniform  stresses  or 
strains  as  well  as  local  transformation  strains. 
The  latter  can  be  generated,  for  example,  by  a 
change  of  temperature  or  inelastic  flow  of  the 
phases,  and  will  be  examined  in  more  detail  in 
the  next  two  sections.  In  the  present  treatment, 
we  assume  that  some  localized  volumes  of  the 
microstructure  undergo  uniform  transforma¬ 
tion  strains  without  reference  to  their  origin. 


1.14.2.1  Governing  Equations 

A  representative  volume  V  of  the  composite 
aggregate  is  selected  such  that  its  overall  re¬ 
sponse  under  uniform  fields  is  identical  to  that 
of  the  composite.  The  local  stresses  and  strains 
within  V  vary  pointwise,  but  will  be  approxi¬ 
mated  by  piecewise  uniform  fields  over  subvo¬ 
lumes  Vn  r  =  1,2,. . .  Q,  such  that  V  =  £F,.,  or 
£cr  =  1,  where  cr  =  Vr/V  is  the  volume  frac¬ 
tion.  The  number  of  subvolumes  Q  indicates 
the  degree  of  refinement  of  the  local  fields.  For 
example,  a  crude  approximation  of  the  local 
stress  and  strain  by  a  piecewise  uniform  distri¬ 
bution  over  two  subvolumes  (one  belongs  to  the 
matrix  and  the  other  to  the  fiber),  provides  the 
averaging  models.  On  the  other  hand,  the  fiber 
and  matrix  can  be  subdivided  into  many  small 
volumes,  which  reside  in  either  phase  as  em¬ 
ployed,  for  example,  in  finite  element  modeling 
of  the  representative  volume  V.  Let  cr,.  and  er 
denote  the  average  stress  and  strain  in  a  sub¬ 
volume  such  that 

V'  =  y  |  ff(x)dK„  er=jF  |E(x)d^  (1) 

Vr  Vr 


where  x  represents  the  coordinates  of  a  point  in 
Vr.  Similarly,  the  overall  uniform  stress  cr  and 
strain  e  can  be  written  as 


<r(x)dV=-J2 


< j{\)i\Vr 


v  vf 

Q  Q 

~  ^  ^  Cr<Tr,  £  —  'y  ^  CrEr 
r- I  r=l 


(2) 
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Consider  elastic  behavior  of  each  homoge¬ 
neous  subvolume  subjected  to  a  uniform  strain 
er,  or  stress  07,  and  simultaneously  undergoing 
uniform  transformation  strain  or  stress  Ar. 
In  the  subsequent  sections,  where  nonlinear 
response  is  considered,  we  will  deal  with  stress 
and  strain  increments  rather  than  total  quanti¬ 
ties.  Consequently,  the  constitutive  equations 
of  the  elastic  phases  or  subvolumes  will  be 
written  in  a  rate  form  as 

o 7  “  L fEf  -t"  Xr,  £r  lVIr(Tr  "t*  |ir  (3) 

where  Lr  and  Mr  =  Lr-1  are  elastic  stiffness 
and  compliance  matrices.  Unless  there  is  a 
change  of  temperature,  these  property  matrices 
remain  constant  within  Vr.  From  Equation  (3), 
the  transformation  stress  and  strain  are  related 
by\r=-L  ri*,r 

Similarly,  the  constitutive  equations  for  the 
composite  aggregate  subjected  to  overall  uni¬ 
form  strain  e,  or  stress  a,  and  undergoing  uni¬ 
form  transformation  strain  |x,  or  stress  X,  are 
written  as 

d  =  Ls  +  X,  8  —  Md  +  (4) 

where  L  and  M  =  LTl  are  overall  elastic  stiff¬ 
ness  and  compliance  matrices,  and  (x  and  X  are 
related  by  X  =  -L|x. 

The  connection  between  the  overall  property 
matrices  and  their  local  counterparts  was  estab¬ 
lished  by  Hill  (1963,  1965b)  in  terms  of  phase 
strain  and  stress  concentration  factors,  denoted 
respectively  by  Ar  and  Br.  In  the  absence  of 
transformation  fields,  the  local  strain  and  stress 
averages  in  the  phases  or  subvolumes  are  writ¬ 
ten  as 

er  =  A  dr  =  Brd  (5) 

If  the  concentration  factors  are  known,  one  can 
utilize  Equations  (2) -(4)  to  find  the  overall 
elastic  stiffness  and  compliance  in  the  form 

Q  Q 

L  =  ^  ^  crLrAr,  M  =  y  ^  crMrBr  (6) 

r=l  r=l 

together  with  the  connections 

£crAr  =  I,  y>B,=I  (7) 

r=l  r=l 

ArM  =  MrBr,  BrL  =  LrAr  (8) 

where  I  is  identity  matrix.  When  Q  =  2,  the 
local  fields  are  determined  as  averages  over 
the  fiber  and  matrix  phases.  In  this  case,  Equa¬ 
tions  (6)  and  (7)  provide  the  concentration 
factors  in  terms  of  the  local  and  overall  moduli 
as 


Ar  =  (L, — L*)  ~ 1  (L — L. v)/ cr 

Br  =  (Mr  -  Ms)  ~ 1  (M  -  M s)/cn  r,s  =f,m  (9) 

In  the  presence  of  local  transformation  fields 
in  the  phases  or  subvolumes,  Vnr  =  1,2,...  Q, 
additional  strains  and  stresses  are  generated 
and  superimposed  on  the  fields  in  Equation 
(5),  to  yield  (Dvorak,  1992), 

Q  Q  . 

er  =  Are  +  y^D^iij,  crr  =  Brcr  +  >  FrA  (10) 

r=  1 

where  and  Frj,  r,  5  =  1,  2,. . .  Q,  are  trans¬ 
formation  influence  functions.  The  kth  column, 
k=  1,  2, . . .  6,  of  Drs  evaluates  the  local  strain  zr 
in  phase  Vn  caused  by  a  uniform  transforma¬ 
tion  strain  component  of  unit  magnitude 
present  in  Vs  under  8  =  0.  Similarly,  the  A;th 
column  of  Frs  evaluates  the  local  stress  ar  in  Vn 
caused  by  a  uniform  transformation  stress  Xjt  of 
unit  magnitude  present  in  Vs  under  a  =  0. 

The  overall  transformation  stress  and  strain 
are  given  in  terms  of  their  local  counterparts  by 
the  generalized  Levin’s  (1967)  formula  (Dvorak 
and  Benveniste,  1992), 

k  =  £cr aJk,  |i  =  cvbJp,  (li) 

r=l  r—  l 

The  mechanical  concentration  factors  and 
the  transformation  influence  functions  depend 
on  the  micromechanical  model  used  in  approx¬ 
imating  the  local  fields.  In  the  sequel,  we  briefly 
describe  two  classes  of  models,  which  can  be 
utilized  to  compute  the  concentration  factors  Ar 
and  Br,  and  the  transformation  influence  func¬ 
tions  and  ¥rs,  r,  s  =  1,  2,...  Q ,  for  the 
subvolumes  Vr.  Under  certain  approximations 
of  the  local  fields,  closed  form  expressions  can 
be  found  for  the  influence  functions  in  terms  of 
the  mechanical  concentration  factors.  Consid¬ 
ering  estimates  derived  with  the  self-consistent 
(Hill,  1965a,  1965b  or  Mori  and  Tanaka,  1973) 
methods,  discussed  in  the  sequel,  Dvorak  and 
Benveniste  (1992)  found  the  following  expres¬ 
sions  for  both  methods, 

D„  =  (I— Ar)(  Lr — L) ~ 1  (§„I  —  csAj )L,  (12) 

Frs  =  (I — B,.)(Mr — M)  “ 1  (8rjI  —  CjBj  )Mj 

r,  s  =  1,2,. . .  Q  (13) 

with  the  connection 

F„  =  Lr(S„I  -  csArB'[  -  Dra)My  (14) 

where  is  the  Kronecker’s  tensor. 
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When  the  local  fields  are  averaged  over  the 
fiber  and  matrix  phases,  0  =  2,  Equations  (12) 
and  (13)  reduce  to  (Dvorak,  1991) 

Dr,  =  (I  —  Ar)(Ly  —  Lm)  “ 1  Ly 

Fw  =  (I-Br)(Mf-Mmrl  M5,  r,s  =  fm  (15) 


1.14.2.2  Micromechanical  Models 

1.14.2.2.1  Two-phase  averaging  models 

These  models  utilize  Eshelby’s  solution 
(Eshelby,  1957)  of  an  ellipsoidal  inhomogeneity 
in  an  infinite  matrix  under  remotely  applied 
uniform  fields  to  estimate  the  average  stresses 
and  strains  in  the  matrix  and  the  fiber  reinfor¬ 
cement.  Hence,  the  number  of  subdivisions  of 
the  representative  volume  V  in  these  models  is 
0  =  2.  Starting  with  the  solution  of  a  transfor¬ 
mation  problem  in  which  a  region  (referred  to 
as  the  inclusion),  in  an  infinite  homogeneous 
isotropic  elastic  medium  undergoes  a  “sponta¬ 
neous  change  of  form,”  or  eigenstrain  p, 
Eshelby  (1957)  established  that  if  the  inclusion 
is  ellipsoidal,  then  the  stress  and  strain  within 
the  inclusion  are  uniform,  but  not  necessarily 
coaxial  with  jjl.  This  important  result  led 
to  the  solution  of  several  related  problems.  In 
the  inhomogeneity  problem,  an  ellipsoid  of  the 
same  shape  and  size  as  the  untransformed  in¬ 
clusion,  but  from  a  different  material,  is  em¬ 
bedded  in  the  matrix,  and  a  uniform  strain  field 
is  applied.  Solution  of  this  problem  is  provided 
by  the  first  inclusion  problem  in  which  an 
“equivalent  inclusion,”  subjected  to  a  certain 
transformation  strain  and  the  given  applied 
uniform  strain  field,  replaces  the  inhomogene¬ 
ity  without  altering  the  stresses  or  strains  any¬ 
where. 

Let  p  be  the  unknown  transformation  strain 
to  be  applied  to  the  equivalent  inclusion,  which 
has  the  same  material  properties  as  the  sur¬ 
rounding  matrix.  The  constrained  uniform 
strain  field  in  the  inclusion  is  given  by  Eshelby 
(1957)  as 

ec  =  Sfju  (16) 

where  S  is  a  constant  matrix  which  depends  on 
the  shape  of  the  ellipsoidal  inclusion  and  Pois¬ 
son’s  ratio  of  the  isotropic  matrix  material.  If  a 
uniform  strain  field  s  is  superimposed,  the 
strain  in  the  inclusion  is  given  by  (s  4*  ec), 
and  the  stress  by  Lj(e  +  ec  -  p),  where  L!  is 
elastic  stiffness  of  the  matrix  and  inclusion.  On 
the  other  hand,  the  stress  in  the  inhomogeneity, 
with  stiffness  L2,  is  L2e2,  where  e2  =  e  4-  ec  is 
the  total  strain.  The  transformation  strain  in 


the  equivalent  inclusion  is  then  found  by  equat¬ 
ing  the  stresses  in  the  inclusion  and  the  inho¬ 
mogeneity,  utilizing  Equation  (16), 

Lt(s  4  Sp-p)  =  L2  (e  4  Sp)  (17) 

and  the  strain  in  the  inhomogeneity  is  found  as 
e2  =  eSp.  Alternately,  e2  can  be  found  directly 
from  Equations  (16)  and  (17)  after  rewriting  the 
latter  as  L{(e2  — p)  =  L2e2.  The  result  is 

e2  =  [I-SLr'(L,  -L2)]-'e  (18) 

This  result  has  been  utilized  in  several  models 
to  determine  the  strain  and  stress  concentration 
factors  in  two-phase  composites.  Considering  a 
dilute  concentration  of  ellipsoids  of  Material  2 
(the  reinforcement)  into  a  matrix  of  Material  1, 
the  strain  concentration  factor  in  the  reinforce¬ 
ment,  A2,  is  given  by  the  coefficient  matrix  in 
Equation  (18)  (Christensen,  1979).  The  corre¬ 
sponding  stress  concentration  factor,  B2,  can 
then  be  found  from  Equation  (8).  Hence, 

A2  =  [I  — SLf 1  (Li-L2)]”1 

B2  —  L2[I  —  SLf^Li  -L2)]-!M  (19) 

where  M  is  the  overall  compliance  of  the  com¬ 
posite  aggregate.  If  the  Eshelby  matrix  S  is 
known,  Equation  (19)  can  be  used  to  find  A2, 
and  Equation  (7),  with  0  =  2,  provides  Ah  the 
matrix  strain  concentration  factor.  The  stress 
concentration  factors  are  then  found  from 
Equations  (19)  and  (7),  where  the  overall  com¬ 
pliance  is  given  by  M  =  L"1,  and  the  overall 
stiffness  L  is  computed  from  Equation  (6)  with 
0  =  2. 

Hill  (1965b)  gives  a  variant  form  of  the  con¬ 
centration  factors  of  Equation  (19)  in  terms  of 
the  stiffness  and  compliance  matrices  of  the 
cavity  containing  the  inhomogeneity.  His  pro¬ 
cedure  leads  to  the  following  form  of  the 
Eshelby  matrix  S 

S  =  PLi  =  I  — M,Q  (20) 

where  P  and  Q  are  related  by  PLj  +  Q  =  I, 
and  depend  on  the  shape  of  the  ellipsoidal 
inclusion  and  elastic  moduli  of  the  matrix.  In 
this  case,  the  concentration  factors  are  written 
as 

Ar  =  [I-P(L,— L,.)]-' 

Br  =  [I-Q(M,-Mr)]-';  r=  1,2  (21) 

We  note  that  for  the  dilute  approximation, 
where  the  reinforcement  causes  only  a  small 
perturbation  of  the  overall  field,  Aj  =  B,  =  I. 
Evaluation  of  the  concentration  factors  and  the 
overall  moduli  is  thus  reduced  to  evaluation  of 
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P  or  Q.  Forms  of  P  for  various  shapes  of 
ellipsoidal  inclusions  can  be  found,  for  exam¬ 
ple,  in  Walpole  (1967,  1969),  Kinoshita  and 
Mura  (1971),  and  Mura  (1982).  Matrix  P  eval¬ 
uated  for  continuous,  parallel  cylindrical  fibers 
embedded  in  a  transversely  isotropic  matrix  is 
given  in  Section  1.14.8. 

Interaction  of  the  inhomogeneities,  which  is 
neglected  in  the  dilute  approximation,  has  been 
considered  in  both  the  self-consistent  (Bu- 
diansky,  1965;  Hill,  1965a,  1965b)  and  the 
Mori-Tanaka  (Mori  and  Tanaka,  1973)  meth¬ 
ods  in  two  different  ways.  In  the  self-consistent 
method,  the  reinforcement  is  embedded  in  a 
homogeneous  medium,  which  has  the  effective, 
but  yet  unknown,  properties  of  the  composite. 
In  this  case,  Equation  (21),  with  Lj  =  L  and 
Mi  =  M  where  L  and  M  are  overall  stiffness 
and  compliance  matrices,  give  the  concentra¬ 
tion  factors; 

Ar  =  [I-P(L  — Lr)]~' 

Br  =  P-Q(M-Mr)]_1,  r  =  1,2  (22) 

Since  the  overall  properties  are  not  known  and 
depend  on  the  concentration  factors,  Equation 
(22)  is  implicit  and  must  be  solved  iteratively  for 
the  concentration  factors,  or  the  overall  mod¬ 
uli. 

The  Mori-Tanaka  method,  on  the  other 
hand,  provides  explicit  forms  for  the  concentra¬ 
tion  factors.  The  solution  is  found  in  two  steps 
(Benveniste,  1987).  First,  the  reinforcement  L2 
is  embedded  in  the  matrix  Lb  and  a  uniform 
displacement  field,  derived  from  the  matrix 
strain,  is  applied  at  the  remote  boundary.  This 
provides  a  partial  strain  concentration  factor, 
T,  which  defines  the  uniform  strain  in  the 
reinforcement  in  terms  of  the  matrix  strain. 
An  equivalent  formulation  can  be  found  by 
subjecting  the  remote  boundary  to  tractions 
derived  from  the  matrix  stress,  and  a  partial 
stress  concentration  factor,  W,  is  found.  Using 


Equation  (21),  we  find 

e2  =  Ts,,  T  =  [I-P(L1-L2)]-1  (23) 

(T2  =  W<fi,  W  =  P-Q(M,-M2)]-'  (24) 

Finally,  Equation  (2)  is  utilized  to  obtain  the 
strain  and  stress  concentration  factors,  for  the 
reinforcement  and  the  matrix,  respectively,  as 

Ai  —  (cj  +  c2T)-1,  A2  =  TA!  (25) 

—  (ctI  +  c2W)-1,  B2  =  WBj  (26) 


Given  the  concentration  factors,  the  overall 
elastic  stiffness  and  compliance  matrices  for 
two-phase  composites  are  determined  from 


Equation  (6)  with  Q  =  2.  Overall  elastic  moduli 
derived  from  the  self-consistent  and  Mori-Ta¬ 
naka  models  for  fibrous  composites  are  given  in 
Section  1.14.8. 


L14.2.2.2  Periodic  array  models 

In  this  class  of  models,  the  actual  material 
microgeometry  of  the  composite  is  replaced 
with  a  certain  periodic  approximation.  Under 
overall  uniform  fields  and  a  uniform  tempera¬ 
ture  change,  the  local  fields  possess  certain 
symmetric  features  such  that  a  unit  cell  can  be 
selected  for  evaluation  of  the  local  fields  and  the 
overall  response  (Iwakuma  and  Nemat-Nasser, 
1983;  Dvorak  and  Teply,  1985;  Teply  and 
Dvorak,  1988;  Nemat-Nasser  and  Hori,  1993; 
Walker  et  al. ,  1994;  Moulinec  and  Suquet,  1994; 
Buryachenko,  1996).  As  an  application  to  fi¬ 
brous  composites,  we  show  in  Figure  1  a  high 
contrast  micrograph  of  the  transverse  plane  of  a 
boron-aluminum  composite  consisting  of  sev¬ 
eral  monolayers,  and  its  idealization  by  a  per¬ 
iodic  hexagonal  array  of  fibers  (Dvorak  and 
Teply,  1985;  Teply  and  Dvorak,  1988).  Several 
other  configurations  of  periodic  dispersions  of 
fibers  in  the  transverse  plane  of  unidirectional 
composites  can  be  selected  (Brockenbrough 
et  al. ,  1991).  We  note  that  the  hexagonal  fiber 
arrangement  provides  rotational  symmetry 
about  the  fiber  longitudinal  axis,  while  idealiza¬ 
tions  with  square  and  rectangular  array  possess 
only  orthotropic  symmetry. 

In  any  case,  a  unit  cell  is  selected  and  sub¬ 
divided  into  small  volumes,  Vr,  r  =  1,  2,. . .  Q 
such  that  each  subvolume  belongs  either  to  the 
matrix  or  to  the  fiber.  Evaluation  of  the  piece- 
wise  uniform  local  stresses  and  strains  in  the 
subvolumes  under  overall  uniform  fields  and 
uniform  temperature  change  can  proceed  in 
two  different  ways.  In  the  first  approach,  the 
overall  stress  and  temperature  time  rates  are 
applied,  and  the  solution  is  found  by  the  finite 
element  method.  The  overall  loading  path  is 
either  known  a  priori ,  or  derived  from  struc¬ 
tural  interaction  of  the  unidirectional  lamina 
under  consideration  with  other  plies  as  encoun¬ 
tered  for  example  in  fibrous  laminates  (Bahei- 
El-Din  et  al ,  1998).  Next,  nodal  forces  that 
equilibrate  the  applied  overall  stresses  are  com¬ 
puted.  Consider,  for  example,  the  unit  cell  de¬ 
rived  from  the  periodic  hexagonal  array  model 
(Dvorak  and  Teply,  1985;  Teply  and  Dvorak, 
1988)  shown  in  Figure  2.  The  displacement 
boundary  conditions  consist  of  the  six  con¬ 
straints  indicated  in  Figure  2  to  eliminate  rigid 
body  motions,  multipoint  constraints  at  the 
triangular  boundary  derived  from  the  assumed 
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(a)  (b) 


Figure  1  Transverse  cross-section  of  a  fibrous  composite:  (a)  high  contrast  micrograph  of  a  B/Al 
composite,  Cf  =  0.45,  and  (b)  idealization  with  the  periodic  hexagonal  array  model. 


periodic  geometry  of  the  microstructure,  and 
generalized  plane  strain  boundary  conditions. 
Assuming  a  linear  displacement  field  in  an 
equivalent  macroscopically  homogeneous 
volume  V,  the  method  of  virtual  work  can  be 
used  to  compute  the  nodal  forces,  pA,  k-  1,2, 

. . 6,  applied  at  the  independent  degrees  of 
freedom  shown  in  Figure  2,  from  the  overall 
stresses.  The  result  is 

h  h 

P\  ~  —-022,  Pl~—  ^023 

h\  h  1 

/>3  =  -yC33+-a23-l(73,  {27) 

h\  h  1 
p*  =  TCT33  +  2^a23  “  i 

Z?5  =  an,  p6  ~  <y\2 

where  £  =  \/3,  and  h  is  the  length  of  the  unit  cell 
in  the  axial  direction  x\,  which  is  to  be  selected 
such  that  the  largest  aspect  ratio  of  the  finite 
elements  is  in  the  order  of  10. 

In  the  second  approach,  the  unit  cell  is 
treated  as  an  aggregate  of  subvolumes  V„ 
r  =  1,  2,. . Q,  and  Equation  (10)  is  solved  for 
the  local  stresses  and  strains.  Since  the  local 
transformation  fields  in  Equation  (10)  may 
depend  on  the  total  or  incremental  stress  or 
strain,  the  form  of  the  governing  equations 
derived  from  Equation  (10)  for  the  local  fields 
vary  according  to  the  underlying  constitutive 
law  of  the  phases.  This  approach,  known  as  the 
transformation  field  analysis,  has  been  applied 
by  Dvorak  et  ai  (1994)  for  viscoelastic,  elastic- 
plastic,  and  viscoplastic  phases,  and  will  be 
described  in  the  subsequent  sections.  In  pre¬ 
paration  for  this  solution  method,  the  concen¬ 


tration  factors  Ar  or  Br,  and  the  transformation 
influence  functions,  Drs  or  F„y,  r,s  =  1,  2,  ...  Q, 
which  depend  on  elastic  properties  of  the 
phases  and  the  assumed  microgeometry,  must 
be  first  evaluated  for  the  selected  subdivision  of 
the  unit  cell. 


V2H; 


Figure  2  Geometry,  constraints,  and  loading  of 
the  unit  cell  derived  from  the  periodic  hexagonal 
array  model. 
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As  suggested  by  Equation  (10),  the  respective 
columns  of  the  concentration  factors  Ar  or  Br, 
represent  the  strains,  or  stresses,  computed  in 
subvolume  Vr  under  unit  overall  strain,  or 
stress,  components,  while  all  subvolumes  are 
free  of  transformation  fields.  Specifically,  the 
fcth  column,  k  =  1,2,. . .,  6,  of  the  stress  concen¬ 
tration  matrix  B„  r  =  1,  2,. . .  Q,  is  given  by  the 
stresses  caused  in  Vr  by  an  overall  stress  com¬ 
ponent  Gjc  —  1.  For  example,  the  first  column  of 
B,  is  given  by  the  stresses  caused  in  Vr  by  the 
overall  stresses  an  =  1,  a,y  =  0,  ij  ^  1.  Substi¬ 
tuting  these  stresses  into  Equation  (27)  provides 
the  nodal  forces 

Pi  =  Pi  -  P3  =  Pa  =  0,  />5=1>  Pe  =  0  (28) 

which  are  applied  to  the  unit  cell  of  Figure  2  to 
compute  the  stress  concentration  factors  using 
the  finite  element  method. 

Similarly,  Equation  (10)  suggests  that  the 
columns  of  the  transformation  matrix,  Dr5, 
r,s  —  1,  2,...,  2,  represent  the  strains  caused 
in  subvolume  Vr  by  unit  transformation  strains 
applied  to  VSi  while  the  unit  cell  is  fully  con¬ 
strained.  Specifically,  if  a  transformation  strain 
component  =  l,  k-  1,  2,. . .,  6,  is  applied  to 
Vs,  the  resulting  strains  in  subvolumes  Vr, 
r=  1,  2,...,  2,  represent  the  kth  column  of 
the  strain  transformation  factor  D^.  Solution 
of  this  problem  can  be  found  by  the  finite 
element  method.  Let  denote  the  transforma¬ 
tion  strain  vector  applied  to  Vs.  Using  the 
principle  of  virtual  work,  and  assuming  linear 
variation  of  the  displacements  within  each  sub¬ 
volume,  the  equivalent  nodal  loads  are  given  by 

p— (29) 

where  L5  is  elastic  stiffness  matrix  of  the  mate¬ 
rial  in  Vs  and  matrix  B  relates  the  nodal  dis¬ 
placements  to  the  element  strains.  The  local 
strains  found  in  the  elements  of  the  unit  cell 
under  the  nodal  forces  given  by  Equation  (29) 
provide  a  column  of  the  strain  transformation 
matrices  Drs ,  which  corresponds  to  the  unit 
transformation  strain  component  found  in  [is . 
A  total  of  62  similar  problems  are  to  be  solved 
to  compute  all  six  columns  of  the  transforma¬ 
tion  influence  functions.  Flowever,  internal 
symmetry  of  the  unit  cell  in  the  transverse 
plane  can  be  utilized  to  simplify  this  process 
(Dvorak  et  al. ,  1994).  Evaluation  of  the  stress 
influence  functions  ¥rs  follows  a  similar  proce¬ 
dure  in  which  a  unit  transformation  stress  com¬ 
ponent  Xk=  1,  k  =  1,  2,. . 6,  is  applied  to  V s. 

Under  isothermal  loading  conditions,  the 
concentration  factors  and  transformation  influ¬ 
ence  functions  remain  constant,  and  are  com¬ 
puted  prior  to  the  analysis  of  the  representative 


volume.  On  the  other  hand,  change  of  the 
elastic  moduli  with  temperature  requires  re- 
evaluation  of  all  the  factors  in  the  course  of 
the  nonlinear  solution.  In  this  case,  efficiency  of 
the  transformation  field  analysis  method  for 
periodic  array  models  may  be  reduced. 


1.14.3  ELASTIC-PLASTIC  ANALYSIS 

We  now  consider  the  solution  of  composite 
aggregates  in  which  behavior  of  the  phases 
under  mechanical  loads  is  nonlinear,  but  rate- 
independent,  beyond  an  elastic  limit  or  yield 
point.  The  models  described  above  for  elastic 
phases  have  been  utilized  in  evaluation  of  the 
local  stresses  and  strains  for  elastic-plastic 
phases.  This  is  described  in  this  section  with 
results,  but  first  a  brief  summary  of  the  elastic- 
plastic  constitutive  equations  of  the  homoge¬ 
neous  phases  is  given. 


1.14.3.1  Plasticity  of  Homogeneous  Materials 

We  consider  a  representative  volume  of  a 
phase  material  of  the  composite  and  regard  it 
as  elastically  homogeneous  and  isotropic  on  a 
macroscale.  As  usually  postulated  for  metals, 
and  verified  by  experiments,  we  assume  the 
existence  of  a  yield  surface  in  the  stress  space, 
which  is  the  locus  of  stress  points  that  can  be 
reached  by  purely  elastic  loading  excursions 
from  the  current  stress  state.  Furthermore, 
load  excursions  which  go  beyond  the  yield 
stress  states  cause  the  yield  surface  to  distort, 
translate,  and  deform  isotropically  such  that 
the  current  stress  point  is  always  contained  by 
the  yield  surface.  Neglecting  distortion  of  the 
yield  surface,  and  assuming  that  its  isotropic 
deformations  are  caused  primarily  by  change  of 
the  yield  stress  with  temperature,  the  current 
yield  surface  can  be  written  in  the  general  form 

/{(a  —  a),7)  =  0  (30) 

where  T  is  the  current  temperature,  and  a  de¬ 
notes  the  position  of  the  center  of  the  yield 
surface.  In  the  absence  of  prior  inelastic  defor¬ 
mation,  a  =  0,  and  Equation  (30)  defines  the 
initial  yield  surface. 

Assuming  that  plastic  deformations  do  not 
affect  the  material  symmetry  found  in  the  elas¬ 
tic  state,  specific  forms  of  the  yield  function  in 
Equation  (30)  are  typically  written  in  terms  of 
stress  invariants.  For  example,  the  Mises  form 
of  the  yield  criterion  is  a  function  of  the  second 
invariant  of  the  stress  deviation  tensor,  a,  and  is 
written  as 
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f=l{i-a):{i-a)-Y2{T)  =  0  (31) 

where  a  is  the  center  of  the  yield  surface  in  the 
deviatoric  stress  space,  and  Y  is  the  yield  stress 
in  simple  tension,  the  magnitude  of  which  may 
depend  on  temperature.  In  Equation  (31)  we 
used  the  notation  (a:b)  to  denote  the  inner 
product  of  second  order  tensors  ay  and  by. 

If  the  strain  decomposition  suggested  by 
Equation  (4)  into  elastic  and  transformation 
parts  is  assumed,  the  stress-strain  relation  is 
written  in  a  rate  form  as 

&  =  L(e  —  kP  —  at),  e  =  Md  +  kP  +  6 it  (32) 

where  a  is  the  coefficient  of  thermal  expansion, 
and  et  is  the  plastic  strain  increment.  The  latter 
is  a  function  of  the  current  stress  state  and  the 
applied  stress  increment.  Whether  or  not  a 
plastic  strain  increment  will  develop  when  a 
stress  increment  is  applied  from  a  current 
state,  which  satisfies  the  yield  function  depends 
on  the  loading  direction.  Under  thermomecha¬ 
nical  load  increments,  the  possible  states  found 
under  loading  from  the  current  yield  surface 
are: 

i'  =  0for/*0'  +  m 
=>  elastic  state  or  unloading 


mercial  aluminum,  and  by  Dvorak  et  ai  (1988) 
on  6061  aluminum,  indicated  that  the  yield  sur¬ 
face  translates  in  the  direction  of  the  applied 
stress  increment  to  accommodate  the  stress 
point.  The  magnitude  of  the  translation,  how¬ 
ever,  must  be  adjusted  if  the  yield  stress  varies 
with  temperature.  Under  thermomechanical 
loading,  the  Phillips  law  is  written  as  (Bahei- 
El-Din,  1990) 

a  =  ya  (37) 

and  the  scalar  y  is  found  by  satisfying  the 
consistency  condition  in  Equation  (36).  For 
the  Mises  yield  surface  in  Equation  (31),  the 
result  is 

nmpt 

where  Y(T)  =  dY/dT. 

The  direction  of  the  plastic  strain  increment 
is  established  with  the  normality  rule,  a  conse¬ 
quence  of  Drucker’s  (1952)  postulate  for  a 
stable  material, 

e'  =  rf  (39) 

The  magnitude  of  the  scalar  f  is  evaluated  from 
Ziegler’s  (1959)  equality, 


e'’  =  0  for  /=  0,  (g):*  +  ^f  =  0 

=>  neutral  loading 

■r  0  for  /«  0,  (g):*  +  |^>0 

plastic  loading 


(34) 


(35) 


(40) 


For  the  Mises  yield  surface  in  Equation  (31), 
Equations  (39)  and  (40)  give 


Stress  states  which  fall  outside  the  elastic 
domain  bounded  by  the  yield  surface  in 
Equations  (30)  or  (31)  are  not  allowed.  Conse¬ 
quently,  the  yield  surface  expands  isotropically 
and/or  translates  to  contain  the  stress  points 
that  satisfy  Equation  (35).  Any  hardening  rule 
postulated  for  the  yield  surface  must  conform 
to  Prager’s  (1955)  consistency  condition 

(K> 

written  here  for  kinematic  hardening,  while 
isotropic  changes  in  the  yield  surface  are  caused 
by  change  of  the  yield  surface  with  temperature. 
Evolution  laws  for  the  center  of  the  yield 
surface  a  were  suggested  early  on  by  Prager 
(1955)  and  modified  by  Ziegler  (1959).  More 
recent  experiments  conducted  by  Phillips  and 
co-workers  (Phillips  et  ai ,  1972;  Phillips  and 
Moon,  1977;  Phillips  and  Lee,  1979)  on  com¬ 


where  H(T)  is  plastic  tangent  modulus  of  the 
stress-plastic  strain  curve,  determined  from  a 
simple  tension  test,  and  defined  for  multidimen¬ 
sional  stress  state  as 


Zp 

tp  =  :  #)* 


(42) 


and  n  is  the  unit  normal  to  the  yield  surface  at 
the  current  stress  point.  For  the  Mises  yield 
surface  in  Equation  (31),  n  is  given  by 


1 


<511 


A  22  0  33  2^23 


A 


=  a  —  a 


2^31 


20  12 


T 


(43) 


The  constitutive  law  described  above  can 
now  be  cast  in  one  of  two  forms  to  be  utilized 
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in  the  micromechanical  models  of  Section 
1.14.2.2.  In  the  first  form,  the  stress-strain  re¬ 
lations  are  written  in  terms  of  instantaneous 
moduli.  Hence,  for  thermomechanical  loading, 
we  write 

&  =  £z  +  £T,s  =  Mg  +  mT  (44) 

where  £  and  M  are  instantaneous,  elastic-plas¬ 
tic  stiffness  and  compliance  matrices,  and  l  and 
m  are  instantaneous  thermal  stress  and  strain 
vectors.  Assuming  additive  decomposition  of 
elastic  and  inelastic  strains,  these  property  ma¬ 
trices  can  be  written  from  the  constitutive  equa¬ 
tions  described  above,  for  a  Mises  yield  surface, 
as  (Bahei-El-Din,  1990) 


t  =  L  +  £p  =  M~\  M  =  M  +  MP  (45) 


(46) 


t  =  1  +  ip  =  —  £  m,  nt  =  a  +  tttP  (47) 


(48) 

where  G  is  elastic  shear  modulus. 

Alternately,  the  stress-strain  equations  can 
be  written  in  the  equivalent  form  suggested  in 
Equation  (4), 


<r  =  Le  4-  X,  e  =  Mct  +  (I 

(49) 

where 

L  =  M”1,  X  =  —  L|i, 

(50) 

|X  =  8  —  se  —  Mp&  +  (a  +  mP)T 

(51) 

For  a  Mises  yield  criterion,  Equation  (51)  pro¬ 
vides 

'fi  H  T)n  +  0lT 

(52) 

1.14.3.2  Plasticity  of  Composite  Materials 

The  micromechanical  models  of  elasticity 
described  in  Section  1.14.2  can  now  be  com¬ 
bined  with  the  phase  constitutive  equations  of 
plasticity  to  compute  the  overall  response  and 
local  stress  and  strain  fields  in  composite  mate¬ 
rials.  The  method  depends  on  the  micromecha¬ 
nical  model  used.  First,  consider  two  phase 


averaging  models,  Section  1.14.2.2.1,  and 
assume  that  both  the  fiber  and  matrix  phases 
are  elastic-plastic  with  instantaneous  properties 
given  by  Equations  (45)-(48).  In  analogy  with 
Equations  (5)-(9),  stress  and  strain  increments 
in  the  phases  can  be  written  for  mechanical 
loading  as 

kr  =  A,z,  ov  =  $r6r,  r=/,m  (53) 

where  Ar  and  %  are  instantaneous  strain  and 
stress  concentration  factors,  given  by 

Ar  =  (£r  -£s)~l(£-  £s)/cr 

%  =  (Mr  ~  Msy " 1  ( M  -  Ms)/cn  r,s  =  f,m  (54) 

The  overall  instantaneous  stiffness  and  compli¬ 
ance  matrices  £  and  M  are  given  by 

£  Cf£fAf  "I-  cm£mAm 

M  —  CfM/Bf  +  cmMm5$m  (55) 

Since  the  response  of  the  material  is  a  func¬ 
tion  of  the  past  loading  history,  the  instanta¬ 
neous  moduli  are  not  known  a  priori ,  and  the 
concentration  factors  cannot  be  evaluated  from 
Equation  (54).  Instead,  averaging  models  de¬ 
rived  in  Section  1.14.2.2.1  for  elastic  phases  are 
used  to  compute  the  concentration  factors, 
after  replacing  the  phase  elastic  properties  by 
the  instantaneous  moduli.  We  recall  that  the 
averaging  models  described  in  Section 
1.14.2.2.1  represent  various  interpretations  of 
the  concentration  factors  given  in  Equation 
(19),  which  was  derived  from  Eshelby’s  equiva¬ 
lent  inclusion  problem.  When  the  matrix  in¬ 
stantaneous  properties  are  used  to  solve  the 
elastic-plastic  problem,  we  are  effectively  sol¬ 
ving  successive  problems  of  the  Eshelby  type 
assuming  that  the  matrix  is  elastically  anisotro¬ 
pic.  In  this  case,  evaluation  of  the  Eshelby 
matrix  S  can  be  performed  only  numerically 
(Ghahremani,  1977). 

The  above  approach  has  been  widely  imple¬ 
mented  for  both  isothermal  and  thermomecha¬ 
nical  loads  (see,  for  example,  Bahei-El-Din, 
1990  and  Lagoudas  et  al ,  1991).  The  method 
however  is  known  to  produce  overall  and  phase 
strain  predictions  that  violate  the  rigorous  Le¬ 
vin’s  formula  in  Equation  (11).  The  source  of 
this  violation  and  related  inaccuracies  of  the 
predictions  are  traced  to  application  of  the 
plasticity  constitutive  law  using  phase  average 
stresses.  A  more  refined  approximation  of  the 
local  fields  is  desired  when  the  constitutive 
behavior  of  the  phases  deviates  from  linearity. 
This  is  offered  by  the  transformation  field  ana¬ 
lysis  (Dvorak  et  ai,  1994). 

The  method  centers  on  Equation  (10)  which 
evaluates  the  local  fields  as  the  sum  of  contri- 
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butions  from  the  overall  fields  applied  to  an 
elastic  aggregate,  and  the  local  transformation 
fields.  If  the  latter  are  caused  by  thermoplastic 
deformation  of  the  phases,  Equation  (51)  is 
applied  to  each  subvolume  to  obtain  the  incre¬ 
ment  of  transformation  strain  corresponding  to 
the  current  stress  state  and  temperature.  Sub¬ 
stituting  Equation  (51)  into  Equation  (10),  we 
find 

Q 

07  =  B r&  -  ^2 Frshs (Mps&s  4  (a,  +rnps)T )  (56) 

5=1 


From  Equation  (56),  the  array  of  local  stress 
increments  {07}  =  {<Ti,  6*2,  •  •  .<tq}  can  be  writ¬ 
ten  as 

{«rr}  =  [diag(I)  +  [F„Mtnr' 

{[Br]a  -  [F„  LJ{«,  +  *?}  t]  (57) 


Equation  (57)  is  written  for  a  Mises  material  as 

i-i-i 


{&r}  = 


diag(I)  4-  j 


Hs 


{ [Br]cr  -  f 


(58) 


Integration  of  Equation  (58)  along  the  loading 
path  (6\7)  provides  the  local  stresses.  The  local 
strain  increments  are  found  from  Equation 
(44),  and  the  overall  inelastic  strain  increments 
from  the  generalized  Levin’s  formula  in  Equa¬ 
tion  (11). 

For  a  Mises  yield  criterion,  the  transforma¬ 
tion  strains  in  the  subvolumes  are  given  by 
Equation  (52).  In  this  case,  Equation  (11)  pro¬ 
vides 


=  e?  4  at  (59) 


the  instantaneous  stress  concentration  factors 
in  sub  volume  Vn  r  =  1,2,. . .  Q ,  are  written  as 

m  =  (diag(I)  +  [F„L ,SM>]\  ~ '  [Br]  (62) 

The  overall  instantaneous  compliance  is  then 
found  as 

Q 

Jt  —  ^  crJir!Br 

o'  <63) 

=  crMr  [diag(I)  +  [F„M^]] [&] 

r=  1 

A  similar  procedure  leads  to 

[A]  =  [diag(I)  +  [D^-err  '[Ar],  (64) 


0 

£  =  L Cr£rA’ 

r—\ 

Q 

=  ^cr£f[dmg(I)  +  [D,,MS£P]]-'[A,] 

r=\ 

(65) 

where  Ar  is  the  instantaneous  strain  concentra¬ 
tion  factor  in  Vn  and  £  is  the  overall  instanta¬ 
neous  stiffness  matrix.  The  overall  thermal 
strain  is  found  from  Equations  (11)  and  (47)  as 

Q 

tn  —  yy  cr  Sf  (qr  4  mf)  (66) 

r=l 


For  a  Mises  yield  criterion,  Equations 
(62)  — (66)  provide 


diag(I)  - 

[At] 


2GS 


[\1  4  HS/3GS 


Dr,M,  [n,nj] 


-l 


(67) 


m  = 


diag(I)  4  2 


n  -i 


H$ 


Fr,L,n,nJ 


[Br]  (68) 


Q 

«  =  ^crBjar  (61) 

r=  I 

where  kp  is  the  overall  plastic  strain  increment, 
and  a  is  the  overall  coefficient  of  thermal  ex¬ 
pansion. 

Equation  (57)  can  be  used  to  derive  overall 
instantaneous  property  matrices  and  thermal 
strains  which  can  be  utilized  in  finite  element 
analysis  of  composite  structures  in  which  the 
elements  are  treated  as  macroscopically  homo¬ 
geneous  with  properties  derived  from  a  micro¬ 
mechanical  model  (Bahei-El-Din  et  a!.,  1981). 
Comparing  Equation  (57)  with  Equation  (53), 


r=l 


2  Gs 


diag(I)  “  {i+h!/2Gs)  [D"M'4nJ]]J  [A, 


M  =  Y^crMr 


dmg(I)+2 


F„L,nvnJ 


[Brl 


(69) 

] 

(70) 

(71) 
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1.14.3*3  Bimodal  Plasticity  of  Fibrous 
Composites 

The  micromechanical  models  discussed 
above  for  elastic-plastic  phases  provide  the 
overall  instantaneous  strains  as  well  as  the 
local  fields,  but  can  also  be  used  to  find  the 
overall  yield  surface.  Consider,  for  example,  a 
two-phase  composite  in  which  the  fiber  is  elastic 
and  the  matrix  is  elastic-plastic  with  a  yield 
criterion  described  by  the  Mises  yield  surface 
(31).  Let  6m  and  am  denote  the  matrix  deviatoric 
stress  and  center  of  the  yield  surface,  respec¬ 
tively.  If  6  and  a  denote  their  overall  counter¬ 
parts,  and  since  both  the  vectors  O  and 
(a—  a)  fall  inside  their  respective  yield  surfaces, 
we  can  write 

iPm  —  am)  =  (72) 

where  Bm  is  the  matrix  stress  concentration 
factor.  Substituting  Equation  (72)  into  Equa¬ 
tion  (31),  after  the  latter  is  rewritten  with  index 
m  to  denote  matrix-related  variables,  the  over¬ 
all  yield  surface  is  found  in  the  deviatoric  stress 
space  as 

f=  l  [Bm(a  -  «)]  :  (Bm(d  -  a)]  -  Y2m(T)=0  (73) 

In  the  overall  stress  space  cr,  the  yield  surface 
in  Equation  (73)  is  given  by 

/=  (a  -  a)T(B^CBm)(a  -  a)  -  Y2m  (T)  =  0  (74) 

where  a  is  the  center  of  the  overall  yield  surface, 
and  C  is  a  symmetric  6x6  matrix  with  the 
following  nonzero  coefficients,  Cn  =  C22  = 
C33  =  1,  C\2  —  C13  =  C23  =  —  1/2,  C44  =  C55 
“  C66  =  3.  Overall  stress  excursions  which  fall 
in  the  stress  domain  enclosed  by  the  ellipsoidal 
yield  surface  in  Equation  (74)  cause  only  elastic 
deformation  of  the  matrix.  Dimensions  of  the 
overall  yield  surface  depend  on  magnitude  of 
the  matrix  yield  stress,  and  the  volume  fraction 
of  the  fiber  and  matrix  as  well  as  their  elastic 
moduli. 

Motivated  by  experimental  results,  the  bimo¬ 
dal  plasticity  theory  (Dvorak  and  Bahei- 
El-Din,  1987)  postulates  that  yielding  of  the 
matrix  according  to  the  criterion  given  by 
Equation  (73)  takes  place  in  fibrous  systems 
only  under  loading  conditions  which  allow 
both  the  fiber  and  matrix  to  participate  in 
carrying  the  applied  load.  Hence,  this  deforma¬ 
tion  mode  is  referred  to  as  the  fiber-dominated 
mode  (FDM).  On  the  other  hand,  the  matrix- 
dominated  mode  (MDM)  assumes  that  the  ap¬ 
plied  load  is  carried  by  the  matrix,  while  the 
fiber  constrains  the  matrix  plastic  deformation 


to  simple  shear  straining  on  planes  that  are 
parallel  to  the  fiber  longitudinal  axis.  The 
matrix-dominated  mode  is  thus  represented  by 
a  variant  of  the  continuum  slip-model. 

The  initial  yield  condition  on  any  potential 
slip  plane  (k)  is  taken  as, 

/(tW)  =  (maxx^j  -x20  =  0  (75) 

where  xns  denotes  the  resolved  shear  stress,  and 
xQ  is  the  matrix  yield  stress  in  simple  shear.  The 
active  slip  system  is  defined  by  the  normal  to 
the  slip  plane,  and  by  the  slip  direction  sj9  so 
that  the  resolved  shear  stress  is 

ro 

Considering  plane  stresses  in  the  x^-plane 
(Figure  3),  where  the  xraxis  is  parallel  to  the 
fiber  longitudinal  axis,  the  components  T!  and 
t2  of  the  resolved  shear  stress  xns  are  given  by 

Ti  =  a2icosP,  x2  =  ^cr22sin  2(3  (77) 


The  maximum  resolved  shear  stress  is  evalu¬ 
ated  from 


For  the  slip  system  on  the  plane  k  =  1,  Equa¬ 
tions  (76)  — (78)  provide 

=^cos-lq2  for  |q|<l,  Pi  =  0  for  |q|^l  (79) 


*3 


Figure  3  Geometry  of  the  conjugate  slip  systems 
of  the  matrix-dominated  deformation  mode. 
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Figure  4  Transverse  cross-section  of  the  initial  yield  surface  of  the  matrix-dominated  deformation  mode. 


where  q  =  cr2l/G22  if  cr22  #  0,  and  the  angle  0t 
between  the  slip  direction  and  the  xraxis  is 
given  by 


tan  0i  =  -sin  p<  (80) 

The  conjugate  system  on  the  plane  k  —  2  in 
Figure  3  is  specified  by  the  angles  p2  and  02  that 
are  related  to  Pi  and  0!  by 


p2  =  n  -  Pi, 02  =  0i  for  |qi  ^  1 

P2  -  O,02  =  0,  for  |q|^l  (81) 


where  0  ^  Pi  ^  n/4  and  0  <  0j  ^  2tt. 

Substituting  the  slip  system  parameters  that 
assure  maxima  of  the  resolved  shear  stress 
under  the  applied  overall  plane  stress  state, 
and  assuming  kinematic  hardening  for  the  ma¬ 
trix,  one  finds  the  overall  MDM  yield  condition 
as 

f  _  (a 21  “021 V 

T0  ) 

V 

fb  =  ^  - 1  =  0  for  |q|  >  1 


^  — 1—0  for  |q|  ^  1  (82) 


These  relations  suggest  that  the  MDM  yield 
surface  in  the  overall  plane  stress  space  is  an 
open  cylinder  with  oval  cross-section  in  the 
^22^2rPlane,  and  generators  parallel  to  the 


<7n-axis.  The  cross-section  of  the  MDM  yield 
cylinder  for  the  initial  state  au  —  a22  ~  0,  is 
shown  in  Figure  4.  Note  that  this  surface  is 
independent  of  phase  moduli  and  volume  frac¬ 
tions. 

The  overall  yield  surface  is  given  by  the 
internal  envelope  of  the  FDM  yield  surface  in 
Equation  (74),  projected  into  the  On,  a22,  <J2\ 
stress  space,  and  the  MDM  yield  surface 
branch  in  Equation  (82).  The  experimental  re¬ 
sults  that  follow  show  that  the  shape  and  posi¬ 
tion  of  the  observed  yield  surfaces  are  closely 
predicted  by  the  bimodal  theory  combined  with 
the  Phillips  kinematic  hardening  rule  (Equa¬ 
tions  (37)  and  (38)). 


1.14.3.4  Experimental  Results  and  Predictions 

In  this  section,  we  present  a  sample  of  the 
experimental  results  obtained  in  an  extensive 
study  of  the  behavior  of  boron-aluminum  fi¬ 
brous  systems  (Dvorak  et  at.,  1988;  Nigam 
et  ai ,  1994a,  1994b),  together  with  predictions 
using  the  methods  outlines  in  Sections  1.14.3.2 
and  1.14.3.3.  The  experiments  were  conducted 
on  thin-walled  tube  specimens  fabricated  by 
diffusion  bonding  of  seven  monolayer  6061- 
Al/B  sheets.  The  fibers  were  aligned  parallel 
to  the  meridians  of  the  tube,  at  a  volume  frac¬ 
tion  Cf  =  0.45.  Figure  5  shows  the  geometry  and 
instrumentation  of  the  tube.  A  cross-section  of 
the  tube  is  shown  in  Figure  1(a). 
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Figure  5  Dimensions  and  instrumentation  of  the  composite  tube  specimen. 


The  specimens  were  loaded  by  axial  force, 
internal  pressure,  and  torsion,  applied  accord¬ 
ing  to  various  routines.  Here  we  consider  a 
specific  loading  regime  in  which  the  specimen 
was  subjected  to  internal  pressure  and  torque. 
This  produced  transverse  normal  stress  ct22, 
longitudinal  shear  stress  a2i  and  axial  normal 
stress  <j\\  =  ±a22.  The  axial  stress  was  caused  by 
the  internal  pressure,  and  was  not  compensated 
for  in  the  experiment.  Consequently,  this  load¬ 
ing  regime  can  be  represented  in  a  stress  plane 
consisting  of  the  shear  stress  a2i  and  the  resul¬ 
tant  normal  stress  yj 4-  a^2  =  ^cr22.  This  ls 
shown  in  Figure  6,  where  the  loading  path  is 
indicated  by  the  diamond  symbols,  labeled 
0,  1 . 11- 

The  overall  initial  yield  surface  and  three 
subsequent  surfaces  that  were  detected  in  the 
experiment  are  shown  in  Figure  6.  The  loading 
path  followed  before  evaluation  of  the  subse¬ 
quent  yield  surfaces  is  indicated  in  the  Figure. 
At  each  point  on  the  yield  surface,  the  yield 
stress  was  defined  by  back  extrapolation  from  a 
small  excursion  into  the  plastic  region,  at  zero 
offset  strain.  The  experimental  yield  surfaces 
were  constructed  by  fitting  the  appropriate 
sections  of  the  bimodal  yield  surfaces,  Equa¬ 
tions  (74)  and  (82).  Apart  from  the  change  in 
the  size  of  the  yield  surfaces,  the  bimodal  theory 


predicts  quite  closely  the  shape  of  the  overall 
yield  surface  in  the  loading  plane  of  this  experi¬ 
ment.  Similar  agreement  between  the  yield  sur¬ 
faces  detected  experimentally  and  those 
evaluated  with  the  bimodal  theory  was  found 
under  several  other  loading  regimes  (Nigam 
et  al. ,  1994a,  1994b). 

Figures  7-9  show  the  comparisons  between 
the  predictions  of  plastic  strains  along  the  path 
of  Figure  6.  The  predictions  were  obtained 
from  the  periodic  hexagonal  array  (PHA) 
model  with  two  different  meshes,  the  matrix 
mode  of  the  bimodal  theory,  and  the  Mori- 
Tanaka  model.  The  matrix  constitutive  equa¬ 
tions  described  in  Section  1.14.3.1  were  used  in 
all  models.  The  comparisons  indicate  that  both 
the  Mori-Tanaka  and  bimodal  theories,  which 
compute  the  plastic  strain  using  the  matrix 
stress  average,  fail  to  approximate  the  actual 
plastic  strains.  In  contrast,  the  PHA  model 
provides  a  reasonable  approximation,  even 
with  a  small  number  of  element  subdivisions. 


1.14.4  VISCOPLASTIC  ANALYSIS 

In  this  section,  we  describe  solution  of  the 
composite  aggregate  when  the  rate  of  loading 
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Figure  6  Loading  path  applied  to  the  B/Al  composite  tube  in  the  transverse  tension-longitudinal  shear 
stress  space,  with  experimental  yield  points  and  matrix-dominated  yield  surfaces. 


Plastic  Normal  Strain,  e22  (1 0  ~2) 

Figure  7  Comparison  of  measured  and  predicted  transverse  normal  stress-plastic  strain  response  of  the 
composite  tube  during  loading  along  the  path  0-11  of  Figure  6. 
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Plastic  Shear  Strain,  28^(1  O'2) 


Figure  8  Comparison  of  measured  and  predicted  longitudinal  shear  stress-plastic  strain  response  of  the 
composite  tube  during  loading  along  the  path  0-1 1  of  Figure  6. 


Plastic  Normal  Strain,  e^lO'2) 

Figure  9  Comparison  of  measured  and  predicted  plastic  normal  strain-shear  strain  path  caused  by  loading 

along  the  path  0-1 1  of  Figure  6. 
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has  a  significant  effect  on  the  nonlinear  re¬ 
sponse  of  the  phases.  We  begin  with  a  brief 
summary  of  the  rate-dependent  constitutive 
equations  for  homogeneous  phases,  and  follow 
with  the  solution  for  composite  materials  sub¬ 
jected  to  thermomechanical  loads. 


1.14.4.1  Viscoplasticity  of  Homogeneous 
Materials 

Several  theories  have  been  postulated  to 
model  the  viscoplastic  behavior  of  homoge¬ 
neous  materials  (Bodner  and  Partom,  1975; 
Liu  and  Krempl,  1979;  Eisenberg  and  Yen, 
1981;  Walker,  1981;  Chaboche,  1989;  Neu, 
1993).  In  a  high  temperature  environment,  cer¬ 
tain  features  of  the  response,  such  as  thermal 
hardening  and  recovery,  must  be  included.  To 
represent  in  situ  behavior  of  the  homogeneous 
phases  of  composite  materials,  the  constitutive 
equations  must  adequately  model  the  response 
under  nonproportional  loads.  Several  of  these 
features  are  captured  by  the  theories  found  in 
the  literature.  Here,  we  briefly  describe  a  rate- 
dependent  theory  with  a  yield  surface  which 
was  developed  by  the  authors  (Bahei-El-Din 
et  al.,  1991),  and  implemented  in  micromecha¬ 
nical  models  of  composite  materials  (Bahei-El- 
Din  et  al .,  1991;  Bahei-El-Din,  1994,  1996) 
subjected  to  thermomechanical  loads.  For  iso¬ 
thermal,  proportional  loading  regimes,  the  the¬ 
ory  degenerates  to  the  equations  proposed  by 
Eisenberg  and  Yen  (1981). 

Assuming  that  the  elastic  response  is  rate- 
insensitive,  a  yield  surface  which  bounds  all 
stress  states  that  produce  pure  elastic  strains  is 
postulated  in  the  general  form  in  Equation  (30), 
or  the  Mises  form  in  Equation  (31).  The  latter  is 
rewritten  here  in  terms  of  the  deviatoric  equili¬ 
brium  stress,  which  represents  the  long-term 
stress  state  that  will  be  attained  by  the  material 
under  constant  strain  induced  by  loading  at  a 
certain  time  rate.  Accordingly,  the  correspond¬ 
ing  yield  surface  is  referred  to  as  the  equilibrium 
surface.  If  both  kinematic  and  isotropic  hard¬ 
ening  are  considered,  the  Mises  equilibrium 
surface  is  written  as 

f=  ~  (a*  -  a)  :  (a*  -  a)  -  ( Y(T)  +  Cl)2=  0  (83) 

where  a  is  the  deviatoric  equilibrium  stress 
tensor,  a  denotes  the  center  of  the  equilibrium 
yield  surface,  Y  ~  Y(T)  is  the  initial  yield  stress 
in  tension,  and  Q  =  Q  (7)  is  an  isotropic  hard¬ 
ening  function. 

For  a  given  deviatoric  stress  tensor  a,  which 
lies  outside  the  yield  surface  in  Equation  (83), 
there  exists  an  equilibrium  stress  a*  which 


satisfies  Equation  (83)  and  is  collinear  with  a 
and  a.  Hence, 


2(y+n)2 


3(a*  —  a)  :  (a*  —  a) 


(a*  —  a)  +  a 


(84) 


The  effective  overstress  A  is  a  measure  of  the 
distance  between  the  actual  stress  point  a  and 
the  equilibrium  point  a*.  It  vanishes  if  the  stress 
point  lies  on,  or  falls  inside  the  yield  surface.  In 
particular, 

A  =  ^  (a  ~  **)  :  (*  ~  if /{a  -  a)  >  0  (85) 


A  =  0  if/(a-a)^0  (86) 

If  the  overstress  is  nonzero,  Equation  (85), 
inelastic  strains  develop  at  the  rate 

®i,1= (vf*(]r>/V’(T))n  (87) 

where  the  functions  k(T)  and  p(T)  are  material 
parameters  and  n  is  a  unit  normal  to  the  yield 
surface  in  Equation  (83)  at  the  current  equili¬ 
brium  stress  point.  The  latter  is  found  from 
Equation  (83)  as  (see  also  Equation  (43)) 

1 

y/^RY+Q.) 

[»ll  <>22  »33  2»23  2<>3|  2<»  12  ]T 

a  —  a*  —  a 


The  time  rate  of  D  is  dependent  on  inelastic 
deformation  and  thermal  recovery  of  the  equi¬ 
librium  yield  surface.  It  is  given  by 

n  =  (o  (T) 

[Qa(T)  ~  ft]£in  -  br(T) |Q  -  Qr(T)\MT)~{)  (89) 

[«-a  (T)\ 


The  functions  a>(7),  Qa(T ),  br(T ),  Qf(7)  and 
p r{T)  are  material  parameters,  and  !in  is  the 
effective  inelastic  strain  rate; 


2=k(T)A^\  =  0 


(90) 


Total  (Qr  =  0),  or  partial  (Qr  #  0)  thermal  re¬ 
covery  is  represented  by  the  second  term  in 
Equation  (89). 

In  analogy  with  Equation  (89),  and  permit¬ 
ting  complete  thermal  recovery  of  kinematic 
hardening,  the  evolution  equation  for  the  center 
of  the  equilibrium  yield  surface  can  be  written 
as 
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<>  =  yv  —  r\r(T) a^Xr{T)  !)d,  6  ~  (<*:a)I/2  (91) 

where  r)r  ( T)  and  %r  (7)  are  material  parameters. 
The  unit  tensor  v  defines  the  direction  of  trans¬ 
lation  of  the  yield  surface  in  the  deviatoric  stress 
space,  and  can  be  specified  according  to  the 
hardening  rules  applied  in  rate-independent 
plasticity  theories.  If  the  Phillips  hardening 


rule  is  selected,  then 

o 

Ik 

II 

(92) 

v  =  n  if  a  =  0 

(93) 

The  factor  y  in  Equation  (91)  is  found  from 
Prager’s  consistency  condition  /  =  0,  when 
translation  of  the  yield  surface  is  specified  by 
the  first  term  in  Equation  (91).  The  result  is 

y  = 

y/2/3k(T)A fW[H(T)  -  <*(T)(na(T)  -  Q)] /(n  :  v) 

(94) 


In  analogy  with  the  equilibrium  yield  surface, 
thermal  recovery  of  isotropic,  as  well  as  kine¬ 
matic,  hardening  of  the  bounding  surface  can 
be  included  in  the  model.  This  is  omitted  here 
for  brevity.  We  only  mention  that  the  recovery 
terms  for  isotropic  and  kinematic  hardening  of 
the  bounding  surface  assume  a  form  similar  to 
those  suggested  above  for  the  equilibrium  sur¬ 
face,  but  with  new  material  parameters. 

With  the  above  constitutive  law,  we  can  now 
write  the  stress-strain  rate  equations  in  the 
form  suggested  in  Equation  (4); 

a  =  Le  +  X,  s  =  Mb-  +  |i  (97) 

where 

L  =  M~\  X  —  — L|jl  (98) 

(X  =  ein  +  at  =  n  +  at  (99) 

and  n  is  given  by  Equation  (88). 


A  two-surface  plasticity  theory  (Dafalias  and 
Popov,  1976)  can  be  used  to  describe  evolution 
of  the  instantaneous  tangent  modulus  H .  In 
analogy  with  Equation  (83),  a  bounding  surface 
is  postulated  in  the  form 

g  =  |(S-«):(S-«)-(?(7)  +  ft)2=0  (95) 

where  a  is  the  deviatoric  bounding  stress  tensor, 
determined  from  equality  of  the  normal  to  the 
equilibrium  yield  surface  n(a*  —  a),  Equation 
(88),  and  the  normal  to  the  bounding  surface 
n(a  —  a),_and  a  is  the  center  of  the  bounding 
surface.  Y(T)  denotes  the  bounding  stress,  given 
by  the  intersection  of  the  asymptotic  part  of  the 
uniaxial  stress-strain  curve  and  the  stress  axis, 
and  Q  —  Q(7)  is  an  isotropic  hardening  func¬ 
tion.  The  instantaneous  tangent  modulus  is 
then  given  by 

H(T)  =  H0(T)  +  h(T) 

/3  y/2 

5=y-0:(a--0) 

where  80  is  the  distance  between  the  equilibrium 
yield  surface  and  the  bounding  surface  at  the 
onset  of  inelastic  deformation.  When  the  equi¬ 
librium  stress  point  lies  on  the  bounding  sur¬ 
face,  5  =  0,  and  the  plastic  tangent  modulus 
77(7)  assumes  the  asymptotic  value  H0(T ), 
which  together  with  the  parameters  h(T)  and 
^(7)  need  to  be  determined  experimentally. 


1.14.4.2  Viscoplasticity  of  Composite 
Materials 

Equations  (97)-(99)  can  be  applied  to  each 
viscoplastic  phase  or  subvolume  of  a  composite 
representative  volume  V  to  compute  the  local 
fields  and  the  overall  response.  This  is  achieved 
by  the  transformation  field  analysis,  which  em¬ 
ploys  Equation  (10)  to  compute  increments  of 
the  local  stresses  as  the  sum  of  contributions 
from  the  overall  load  and  the  local  transforma¬ 
tion  strains.  Substituting  Equation  (99)  into 
Equation  (10),  we  obtain  the  stress  rate  in 
volume  Vr  as 

<rr  =  Br&  -  J2  ( JL(7Xs(T)V„L,ns 

y  ’  (ioo) 

-  T,r=l,2,...Q 

s—  1 

Integration  of  Equation  (100)  along  the  loading 
path  (dr,  0)  provides  the  local  stresses.  The  local 
strain  increments  are  found  from  Equations 
(97)  and  (99),  and  the  overall  inelastic  strain 
increments  from  the  generalized  Levin’s  for¬ 
mula  in  Equation  (11). 

This  approach  is  applicable  to  both  the  aver¬ 
aging  models  ( Q  =  2)  in  Section  1.14.2.2.1,  and 
the  periodic  array  models  ( Q  >  2)  in  Section 
1.14.2.2.2.  Selection  of  a  particular  model  is 
reflected  in  the  form  of  the  stress  concentration 
factors  Br  and  influence  functions  ¥rs.  Dvorak 
et  ai  (1994)  illustrate  the  application  of  this 
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Figure  10  Axial  tension-transverse  tension  path  applied  to  a  SCS6/Ti-15-3  composite. 


method  to  unidirectional  composites  with  a 
viscoplastic  matrix. 

Solution  of  thermo-viscoplastic  problems  of 
periodic  array  models  can  be  also  obtained  by 
the  Finite  element  method.  This  is  described  by 
Bahei-El-Din  (1996)  and  implemented  in  the 
Viscopac  code  (Bahei-El-Din,  1994). 


1.14.4.3  Application 

To  illustrate  the  procedure  described  above, 
we  examine  the  time-dependent  response  of  a 
high-temperature  metal  matrix  composite.  The 
composite  material  consists  of  a  titanium  ma¬ 
trix  (Ti-15-3)  reinforced  with  a  silicon  carbide 
fiber  (SCS6)  at  30  volume  percent,  and  was 
examined  at  482°  C.  At  this  temperature,  the 
fiber  is  assumed  to  be  elastic,  and  the  matrix 
elastic-viscoplastic  with  constitutive  behavior 
described  by  the  constitutive  law  given  in  Sec¬ 
tion  1.14.4.1.  Johnson  et  ai  (1993)  give  material 
properties  of  the  constituents. 

The  unit  cell  derived  from  the  periodic  hex¬ 
agonal  array  model  (Figures  1(b)  and  2),  was 
subdivided  into  16  elements,  10  in  the  matrix 
and  6  in  the  fiber.  The  finite  element  method 
was  then  used  to  compute  the  mechanical  and 
transformation  concentration  factors  Br  and 
Frs.  Equation  (100)  and  the  evolution  equations 
for  the  viscoplastic  matrix  in  Section  1.14.4.1, 
were  integrated  for  a  specified  overall  stress 
path  (Dvorak  et  ai ,  1994). 

The  loading  consisted  of  axial  and  transverse 
normal  stresses,  au  and  cr22.  which  were  com¬ 
bined  as  shown  in  Figure  10  and  applied  at  the 
rates  given  in  Figure  11.  The  corresponding 
overall  strains,  En  and  e22,  are  plotted  in 


Figure  12.  The  effect  of  the  time-dependent 
deformation  of  the  matrix  is  evident,  particu¬ 
larly  at  the  sustained  stress  in  the  time  segment 
4-5.  It  is  seen  that  overall  creep  strains  are 
developed  in  both  the  axial  and  transverse 
directions,  albeit  at  different  rates. 


1.14.5  INELASTIC  LAMINATES 

In  this  section,  we  describe  analysis  of  fibrous 
composite  laminates  with  inelastic  phases  using 
micromechanical  models  for  the  unidirectional 
plies.  The  approach  expands  the  transforma¬ 
tion  field  analysis  to  laminates  and  evaluates 
the  instantaneous  response  of  the  laminate  as 
well  as  the  stresses  and  strains  in  the  individual 
plies  and  their  phases.  This  is  achieved  by  re¬ 
garding  the  laminate  as  an  aggregate  of  homo¬ 
geneous  phases  with  constraints  derived  from  a 
micromechanical  model  of  the  ply,  and  the  in¬ 
plane  equi-strain  condition  found  in  laminates. 

First,  we  find  the  lamina  overall  stresses  as  a 
function  of  the  applied  load  and  the  lamina 
transformation  strains.  Next,  we  proceed  to 
evaluate  the  local  stresses  in  the  phases  of  the 
individual  plies  using  two  approaches.  In  one 
approach,  the  transformation  field  analysis  is 
used  to  formulate  rate  equations  for  the  phase 
stresses  in  all  plies  in  terms  of  phase  and  ply 
concentration  factors  and  influence  functions. 
Either  averaging  models  or  periodic  array  mod¬ 
els  can  be  utilized.  In  the  other  approach,  de¬ 
veloped  for  laminates  in  which  the  individual 
plies  possess  a  periodic  microstructure,  the  la¬ 
mina  stresses  arc  applied  to  a  unit  cell  of  the  ply, 
and  the  local  stresses  are  computed  with  the 
finite  element  method. 
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Figure  12  Overall  strain  histories  computed  for  the  loading  path  of  Figures  10  and  11. 


1.14.5.1  Lamina  Stresses 

The  behavior  of  a  symmetric  laminated  plate 
consisting  of  IN  fully  bonded  thin  elastic  plies 
under  thermomechanical  loads  is  considered 
(Figure  13).  Referred  to  a  Cartesian  coordinate 
system,  xk,k  =  1,2,  3,  in  which  the  x\  x2-plane 
coincides  with  the  midplane  of  the  laminate,  in¬ 
plane  membrane  forces  and  the  corresponding 
uniform  stresses,  an,  a22,  and  a12  are  applied, 
together  with  a  uniform  normal  stress,  a33,  in 
the  thickness  direction  x3.  The  out-of-plane 
normal  stress  can  be  found  in  applications 


that  involve  pressure  loading,  such  as  in  fabri¬ 
cation  by  isostatic  pressing,  tr  =  [an,  a22,  ai2] 
lists  the  in-plane  stresses  applied  to  the  lami¬ 
nate,  and  s  =  [sn,  e22>  2ei2]  lists  the  corre¬ 
sponding  laminate  strains.  The  latter  are 
caused  by  the  applied  stresses,  a  and  a33  in 
addition  to  the  transformation  strains  gener¬ 
ated  in  the  individual  plies,  such  as  thermal 
and  inelastic  strains,  which  are  not  recovered 
by  removal  of  the  mechanical  load. 

Assuming  additive  decomposition  of  the  var¬ 
ious  effects  (Dvorak,  1991),  and  adopting  the 
notation  of  Bahei-El-Din  (1992),  the  time  rates 
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Figure  13  Geometry  of  a  symmetric  fibrous  laminate. 


of  the  laminate  in-plane  stresses  and  strains  are 
written  as 

&  =  Ls  4-  ka33  +  X  (101) 

e  =  Md  +  fcr33  +  |ul  (102) 

where  X  —  [Xj  i ,  X,22>  X12]  and  |x  =  [p]  1 ,  jH22->  2^12] 
are  the  in-plane  transformation  stress  and 
strain,  respectively,  and  L  and  M  are  the  elastic 
stiffness  and  compliance  matrices  for  in-plane 
loading.  The  k  vector  lists  the  in-plane  stresses 
caused  in  the  laminate  by  a  unit  out-of-plane 
normal  stress,  cr33  =  1,  in  the  absence  of  both 
the  total  in-plane  strain  £,  and  transformation 
stress  X,  whereas  f  is  the  elastic  compliance 
vector  associated  with  the  out-of-plane  normal 
stress,  a33.  On  the  other  hand,  the  transforma¬ 
tion  strain  |x  represents  the  total  strain  that 
remains  in  the  laminate  after  complete  unload¬ 
ing  to  zero  stress,  and  the  transformation  stress 
X  is  seen  to  represent  the  total  stress  caused  in  a 
fully  constrained  laminate  by  the  transforma¬ 
tion  strain  fx.  Equations  (101)  and  (102)  pro¬ 
vide  the  relations, 

L  =  M-1,  k  =  —Lf,  X=-Ljx  (103) 


and  strain  caused  by  a  uniUni  t-of-plane  normal 
stress,  and  X,  and  jT,  =  —  M/  X,  are  the  ply  in¬ 
plane  transformation  stress  and  strain.  For  a 
transversely  isotropic  ply  with  overall  elastic 
longitudinal,  and  transverse  moduli,  EL  and 
Et,  Poisson’s  ratios,  vL  and  vT,  and  longitudi¬ 
nal  shear  modulus  GL,  the  matrices  found  in 
Equations  (104)  and  (105)  are  given  by  (Bahei- 
El-Din,  1992;  Dvorak  and  Bahei-El-Din,  1995), 


L  ,= 


1 


k  +  m 


Eik  -f  mn 
2  mi 
0 


2nd 
4  km 
0 


0 

0 

p(k  +  m) 


-1 


=  M, 


ki~k  +  m 


r  i 

k  -  m 
0 


(106) 


M,= 

f/  = 


1/  Ei 
-Vl/Ei 
0 

-vl/El 

— vx/£*t 

0 


-vl/Fl  0 

\/Et  0 

0  1  /CL 

=  -M,k/ 


(107) 


In  analogy  with  Equations  (101)  and  (102), 
the  uniform  in-plane  stress  and  strain  rates  of  a 
ply  (/),  i=  1,  2,...,  N,  in  the  local  coordinate 
system  xk,  k  =  1,2,3  (Figure  13),  can  be  written 
as, 


where  k,  I ,  m ,  n ,  and  p  are  Hill’s  (1964)  moduli 
(see  Section  1.14.8),  and  EL~  n  —  fjk,  vL  = 
l/2k,  m  =  Et/2(\  +  vT),CL=Jp. 

When  expressed  in  the  overall  coordinate 
system  xk,  k  =  1,  2,  3,  Equations  (104)  and 
(105)  are  written  as  (Bahei-El-Din,  1992) 


07  =  L /!/  +  k/d33  +  X/  (1 04) 

S/  =  M/d-,-  -f  f/d^  +  pi/  (105) 

where  L,-  and  M,  are  the  elastic  stiffness  and 
compliance  matrices  for  in-plane  loading,  k, 
and  f =  —  M,  k,  are  the  ply  in-plane  stress 


07  =  L/C/  +  k,a33  +  X, 

(108) 

e,  =  Mi&i  +  f/d33  +  pi/ 

(109) 

where 

o-/  =  R,d-„  d33  =  0-33,  e,  =  N/E/ 

(110) 
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K  =  n,-'  = 


L,P„  p,  = 

Nip,- 

(in) 

11 

£ 

=  Nfki  = 

— L/f, 

(112) 

COS2  (p; 

sin2  <p,- 

—  |  sin  2  (p; 

sin2  (p,- 

COS2  <p  j 

^sin  2  (p/ 

_  sin  2(p,- 

-  sin  2cpy 

cos  2(p,- 

(113) 


and  is  the  angle  between  the  local  xraxis  and 
the  overall  xi-axis  (Figure  13). 

The  ply  stresses  in  a  laminate  loaded  by  over¬ 
all  in-plane  stress  ct,  out-of-plane  normal  stress 
ct33,  and  ply  transformation  stress  A„  intro¬ 
duced  by  certain  prescribed  in-plane  transfor¬ 
mation  strains  p,  =  -M,X„  can  now  be 
determined  using  the  transformation  field  ana¬ 
lysis  method.  The  laminate  is  regarded  as  elas¬ 
tic,  and  the  ply  in-plane  stress  is  written  as  the 
sum  of  the  overall  stress  and  local  transforma¬ 
tion  stress  contributions  (Dvorak  and  Bahei-El- 
Din,  1995).  On  the  other  hand,  the  normal 
stress  in  the  thickness  direction  equals  the  ap¬ 
plied  stress.  Hence, 

N 

&i  =  H  (T  +  K,(T33  +  ^  /A/.  O33  =  033  j 

i  =  1,2,  ..AT 

We  note  that  the  lamina  out-of-plane  transfor¬ 
mation  stresses  X33,  X.13,  and  X.23  do  not  neces¬ 
sarily  vanish,  but  they  are  not  introduced  in 
Equation  (114)  since  the  in-plane  equi-strain 
condition  imposed  on  the  perfectly  bonded 
plies  can  be  maintained  under  these  transfor¬ 
mation  stresses  without  introducing  additional 
ply  stresses.  The  H,  and  k,  matrices  are  stress 
distribution  factors  for  in-plane  overall  stresses, 
and  out-of-plane  normal  stress,  respectively, 
and  K,y  is  a  transformation  influence  function. 
The  kth  column  of  matrix  K,y  provides  the  in¬ 
plane  stresses,  ctu,  ct22,  and  cti2  caused  in  la¬ 
mina  (0  by  a  unit  transformation  stress  Xk, 
k  =  1,2,3,  applied  to  lamina  (J)  while  the  overall 
stresses  cr  and  <733  are  absent. 

The  distribution  factors  H,  and  k„  and  the 
influence  coefficients  K,y  are  evaluated  by  rea¬ 
lizing  the  in-plane  strain  compatibility  of  the 
perfectly  bonded  plies,  s  =  e„  and  the  force 
equilibrium  condition,  £  c,  07  =  c,  1  =  1,2,.. 
N,  where  c,-  =  tjt  (see  Figure  13)  is  the  volume 
fraction  of  the  ply.  From  these  conditions,  and 
using  Equation  (114),  one  can  establish  that 

H,  =  L,M,  k,  =  L,(f  -  f,)  (115) 

M  =  L“\  L  =  ^c,Lf  (116) 

i=i 


N 

f  =  — Mk,  k  =  ^c,k, 

1=1 


(117) 


Ky  =  5y  I  —  c,  H,  (118) 

where  8 y  is  Kronecker’s  tensor,  I  is  the  identity 
matrix,  and 


N  N 

£C,H,=I,5>k,=0 


1=1 


1=1 


(119) 


1.14.5.2  Transformation  Field  Analysis 
Method 

We  now  consider  evaluation  of  the  phase 
stress  rates  in  each  lamina  using  the  governing 
equations  presented  in  the  preceding  section, 
when  the  plies  are  subjected  to  the  stress  in 
Equation  (114).  Since  we  are  primarily  con¬ 
cerned  with  transformation  strains  of  thermal 
and  inelastic  type,  we  first  substitute  the  lamina 
transformation  stress  \j  in  Equation  (114)  in 
terms  of  the  eigenstrains  as  X,  =  —  L and  use 
Equation  (111)  to  transform  the  strain  p,  into 
the  lamina  local  coordinates.  Next,  the  ply 
overall  stress  in  Equation  (114)  is  transformed 
into  the  local  coordinate  system  of  the  ply  using 
Equation  (110).  Hence,  the  ply  stresses  in  local 
coordinates  is  written  as 


07  —  R,H,dr  -4-  R  k  <733  —  R,  )  Kyl-  .R  fXj 
M 

^33  =  033,  i=\,2,..N 


(120) 


Given  the  lamina  stress  referred  to  the  mate¬ 
rial  axes,  Equation  (120),  the  local  stresses  in 
the  phases  of  the  individual  plies  are  determined 
using  the  procedure  described  in  Section  1.14.2. 
A  representative  volume  Vj  of  a  unidirectional 
composite  ply  (j)J  =  1,2,. .  .,N,  is  selected  and 
divided  into  subvolumes  Fn,  r|  =  1,2,...,  Q, 
such  that  Vj  =  £  V\.  The  number  of  subvo¬ 
lumes  Q  depends  on  the  micromechanical 
model  selected  for  the  lamina.  For  two-phase 
averaging  models,  Section  1.14.2.2.1,  the  num¬ 
ber  of  subvolumes  is  two  representing  the  ma¬ 
trix  and  the  fiber.  On  the  other  hand, 
representative  volumes  derived  from  periodic 
array  models,  Section  1.14.2.2.2,  may  contain 
several  subvolumes  in  each  of  the  fiber  and 
matrix  phases.  From  Equation  (10),  the  local 
stresses  in  the  subvolumes  of  lamina  (/)  are 
written  as  a  (6  x  1)  vector  in  the  form 
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<  =  +  bP^33  -  £  FPnLn*4 

n=i 

p  =  1,2,  ..g,  /=  1,2,  ..N 


(121) 


The  first  two  terms  represent,  respectively,  the 
local  stress  used  by  the  ply  in-plane  stresses  and 
out-of-plane  normal  stress.  The  (6  x  3)  Bp  ma¬ 
trix,  and  the  (6  x  1)  bp  vector  are  stress  con¬ 
centration  factor  given  by  the  respective 
partitions  of  the  regular  (6  x  6)  stress  concen¬ 
tration  factor  Bp  after  the  columns  of  the  latter 
are  reordered  according  to  the  stress  labels  11, 
22,  12,  33,  13,  23.  The  last  term  in  Equation 
(121)  is  the  local  stress  caused  by  the  phase 
transformation  strains. 

The  lamina  transformation  strain  is  given  in 
terms  of  their  local  counterparts  by  the  Levin’s 
formula  in  Equation  (11),  written  here  for  the 
in-plane  components; 


(122) 


From  Equations  (120)  —  (122),  the  local  stres¬ 
ses  in  the  plies  can  be  written  as 

tfp  =  BpR/H/o-  +  (BpR/K/  +  bp)a33  -  ^ 


Tl  =  l 


-BpR.-EwV 

i=  i 

p  =  1,2,  ..Q,i=  1,2  ,..N 


£<  *n  * 
.n=i 


(123) 


The  first  and  second  terms  in  Equation  (123) 
provide,  respectively,  the  local  stress  caused  by 
the  overall  in-plane  stresses  and  out-of-plane 
normal  stress  applied  to  the  laminate,  while  the 
last  two  terms  are  the  contributions  of  the 
subvolume  transformation  strains  in  all  plies 
to  the  subvolume  p  of  lamina  (/).  The  third 
term  provides  the  local  stresses  due  to  local 
transformation  strains  in  lamina  (/).  The  in¬ 
plane  constraint  e  =  e,  imposed  on  the  lamina 
causes  additional  stresses  in  the  subvolumes  of 
the  plies  when  transformation  strains  fiJn  are 
present  in  other  layers  (/).  This  effect  is  given  by 
the  last  term  in  Equation  (123). 

Considering  a  specific  form  for  the  phase 
transformation  strains,  derived  for  example 
from  thermal  and  viscoplastic  strains  in  Equa¬ 
tion  (99),  the  rate  equations  provided  by  Equa¬ 
tion  (123)  for  the  local  stresses  in  all  plies  can  be 
integrated  along  a  specified  loading  path  (<r, 
a33,  T)  applied  to  the  laminate.  If  the  phase 
elastic  moduli  change  with  temperature,  the 
local  stress  concentration  factors  and  influence 
coefficients  in  all  plies  must  be  updated  within 


the  integration  process.  Wafa  (1994)  gives  ex¬ 
amples  of  laminate  analysis  using  this  ap¬ 
proach. 


1.14.5.3  Finite  Element/Transformation 
Analysis  Method 

In  this  method,  inelastic  analysis  of  laminates 
is  performed  considering  a  detailed  representa¬ 
tion  of  the  microstructure  of  each  ply  using  a 
unit  cell  model.  While  the  analysis  on  the  mi¬ 
croscale  is  performed  with  the  finite  element 
method,  the  stress  rates  for  each  ply  are  ob¬ 
tained  from  a  transformation  field  analysis  of 
the  laminated  plate  (Bahei-El-Din  et  al ,  1998). 
In  this  way,  local  phenomena  related  to  the 
plies  are  incorporated  in  the  finite  element  solu¬ 
tion,  while  the  corresponding  overall  ply  defor¬ 
mation  is  accounted  for  in  the  laminate 
analysis. 

Transformation  field  analysis  of  the  lami¬ 
nates  follows  the  procedure  described  in  the 
preceding  section,  leading  to  Equation  (120). 
This  provides  the  stress  rates,  a\\,  (^2,  cr33,  and 
012',  i=  1,2,...,  N,  for  the  plies,  referred  to 
their  local  coordinate  system.  Response  of  the 
unit  cell  representing  each  lamina  under  these 
stresses  is  computed  with  the  finite  element 
method.  For  example,  if  the  unit  cell  derived 
from  the  PHA  model  is  selected  (Figure  2),  the 
rates  of  nodal  forces  equivalent  to  the  lamina 
stress  rates  are  computed  from  Equation  (27) 
and  applied  at  the  degrees  of  freedom  indicated 
in  Figure  2. 

Considering  viscoplastic  phases,  augmenta¬ 
tion  of  the  finite  element  procedure  with  Equa¬ 
tion  (120)  provides  a  system  of  first  order 
differential  equations  (ODE)  in  the  form 
(Bahei-El-Din  et  al.,  1998), 

yM  =5/(^.Vi,y2,...,^)  (124) 

The  unknown  functions  gj,  j  =  1,  2,. . .,  R ,  are 
identified  with  the  laminate  overall  strain  e,  ply 
stresses  tr,  and  a33  in  the  overall  axes,  and 
transformation  strain  |x„  i  =  1,2,.. .,  N ,  phase 
stress  an  and  transformation  strain  |xn,  r\  =  1, 
2,. . .,  Q  in  all  N  plies.  The  unknown  functions 
also  include  the  nodal  displacements  for  the 
unit  cell  of  each  ply  and  any  internal  variables 
required  to  define  the  rate-dependent  deforma¬ 
tion  of  the  phases,  e.g.  the  overstress  in  Equa¬ 
tion  (85).  Assuming  elastic  response  of  the 
phases  in  the  initial  state,  Equation  (124)  can 
be  integrated  over  a  specified  time  period  using 
an  ODE  solver  that  is  appropriate  for  stiff 
differential  equations  which  are  usually  en¬ 
countered  in  viscoplastic  response  modeled 
with  the  power  law  assumed  in  Equation  (87). 
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Figure  14  Loading  history  applied  to  a  high-temperature  metal  matrix  laminate. 


Figure  15  Overall  axial  strain  computed  for  a  high-temperature  metal  matrix  laminate  subjected  to  the  load 

history  of  Figure  14. 


1.14.5.4  Application 

The  methods  described  above  have  been  used 
to  evaluate  the  overall  strains  and  fiber  axial 
stress  in  a  silicon  carbide-titanium,  (0/±45)s 
laminate  caused  by  in-plane  normal  and  shear 
stresses  applied  in  a  proportional  path  at  565  °C 
(Bahei-El-Din  et  al ,  1998).  The  load-time  re¬ 
cord  is  shown  in  Figure  14.  Material  properties 
of  the  elastic  Sigma  fiber  and  the  viscoplastic 
Timetal-21S  matrix  provided  by  Bahei-El-Din 
and  Dvorak  (1997)  are  used,  and  a  fiber  volume 
fraction  of  0.325  is  assumed. 

The  computed  axial  and  shear  strains  are 
plotted  in  Figure  15  and  16,  and  the  axial  fiber 
stress  in  the  0°-ply  is  shown  in  Figure  17.  The 
figures  compare  the  predictions  obtained  with 
the  Mori-Tanaka  averaging  model  and  the 
periodic  hexagonal  array  model.  In  the  latter, 
a  refined  mesh  of  the  unit  cell  with  48  matrix 
elements  and  24  fiber  elements  was  used  for 
each  ply.  It  is  seen  that  averaging  the  local  fields 


over  the  fiber  and  matrix  phases,  as  modeled  by 
the  Mori-Tanaka  scheme,  underestimates  the 
overall  strains  and  the  fiber  stress  in  compar¬ 
ison  with  the  more  refined  representation  of  the 
local  fields  offered  by  the  finite  element  solution 
of  the  PHA  unit  cell.  Since  axial  deformation  of 
the  laminate  is  dominated  by  the  elastic  0°-fiber, 
the  Mori-Tanaka  estimates  of  the  laminate 
maximum  axial  strain  and  fiber  axial  stress  in 
the  0°-ply  are  smaller  than  the  finite  element 
estimates^by  only  10%.  In  contrast,  a  much 
stiffer  shear  response  is  obtained  with  the 
Mori-Tanaka  model  leading  to  an  estimated 
laminate  maximum  shear  strain  that  is  smaller 
than  the  finite  element  value  by  60%. 


1.14.6  CLOSURE 

Although  efforts  to  estimate  the  overall 
elastic  moduli  of  composite  aggregates  from 
local  properties  date  back  to  the  early  1960s, 
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SCS6/Ti,  cf=  0.325 


Figure  16  Overall  shear  strain  computed  for  a  high-temperature  metal  matrix  laminate  subjected  to  the 

load  history  of  Figure  14. 


Fiber  Axial  Stress  in  0°  Ply  (MPa) 


Figure  17  Fiber  axial  stress  computed  in  the  0°-ply  of  a  high-temperature  metal  matrix  laminate  subjected 

to  the  load  history  of  Figure  14. 


development  of  micromechanical  theories  for 
inelastic  composite  materials  has  seen  major 
advances  only  in  recent  years.  These  have 
been  summarized  in  this  chapter  with  emphasis 
on  fibrous  composites  with  phases  that  exhibit 
plastic  and  viscoplastic  deformation.  Capabil¬ 
ities  of  the  analytical  methods  developed  for 
these  composite  systems  in  predicting  their 
overall  response  under  thermomechanical 
loads  have  been  also  illustrated. 

The  methods  discussed  center  on  evaluation 
of  local  stress  and  strain  concentration  factors 
that  reflect  interactions  of  the  phases  and  their 
in  situ  constitutive  behavior.  Two  classes  of 
micromechanical  models  can  be  utilized  in  com¬ 
puting  the  concentration  factors,  averaging 
models,  which  derive  phase  interactions  from 
solution  of  certain  inclusion  problems  without 


specific  reference  to  the  microgeometry,  and 
periodic  array  models,  which  derive  phase  in¬ 
teractions  from  an  idealized  geometry  of  the 
microstructure.  Representation  of  the  local 
fields  in  these  two  models  is  quite  different.  In 
averaging  models,  the  local  fields  are  averaged 
over  each  phase,  and  thus  the  local  phase  prop¬ 
erties  are  assumed  to  be  homogeneous  within 
the  phase.  Although  valid  for  elastic  phases, 
this  assumption  may  lead  to  large  errors  in 
the  predicted  overall  response  when  behavior 
of  the  phases  is  nonlinear.  In  periodic  array 
models,  a  piecewise  uniform  field  is  selected 
with  a  desired  level  of  refinement.  This  is  ex¬ 
pected  to  provide  more  realistic  predictions, 
albeit  at  a  higher  computational  effort. 

In  plasticity,  formulation  of  the  governing 
equations  for  the  local  stresses  may  proceed  in 
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two  different  ways.  In  one  method,  the  formu¬ 
lation  follows  that  of  elastic  inclusion  pro¬ 
blems,  leading  to  averaging  models,  but  with 
instantaneous  phase  properties.  The  second 
method,  known  as  transformation  field  analy¬ 
sis,  finds  the  local  fields  by  superposition  of  the 
fields  caused  by  the  overall  loads,  and  the  local 
transformation  strains.  The  method  is  applic¬ 
able  to  both  averaging  and  periodic  models, 
and  requires  derivation  of  only  elastic  concen¬ 
tration  and  transformation  factors. 

Viscoplastic  analysis  for  either  averaging  or 
periodic  models  with  the  transformation  field 
method  leads  to  closed-form  rate  equations  for 
the  local  discrete  stress  field.  These  equations 
have  been  expanded  to  include  the  stresses  in 
various  plies  of  a  laminated  structure.  The  re¬ 
sulting  equations  reflect  interaction  of  the 
phases  in  the  individual  plies,  and  mutual  con¬ 
straints  caused  by  bonding  the  plies  together. 

While  the  focus  in  the  present  summary  of 
nonlinear  analysis  of  composites  was  on  micro¬ 
mechanical  models,  other  related  techniques  for 
evaluation  of  the  overall  response  can  be  found 
in  the  literature.  For  example,  several  varia¬ 
tional  methods  for  estimating  effective  behavior 
of  nonlinear  composites  with  random  micro¬ 
structures  have  been  developed  in  the  past  15 
years.  In  elasticity,  this  line  of  inquiry  was 
initiated  by  Hashin  and  Shtrikman  (1962, 
1963)  who  used  special  forms  of  classical  varia¬ 
tional  principles  of  elasticity,  together  with  cer¬ 
tain  polarization  fields  in  a  homogeneous 
comparison  medium,  to  obtain  rigorous 
bounds  on  the  effective  elastic  stiffness  tensor 
of  statistically  homogeneous  composites.  Gen¬ 
eralization  of  the  Hashin-Shtrikman  varia¬ 
tional  principles  to  nonlinear  elastic 
composites  was  first  proposed  by  Talbot  and 
Willis  (1985).  New  variational  principles  using 
a  linear  comparison  composite  were  later  devel¬ 
oped  by  Ponte  Castaneda  (1992),  and  by  Su- 
quet  (1983)  for  power-law  composites.  Several 
other  procedures  and  more  recent  develop¬ 
ments  of  variational  methods  have  been  re¬ 
viewed  and  summarized  by  Suquet  (1997)  and 
Ponte  Castaneda  and  Suquet  (1998). 

Although  they  offer  rigorous  estimates  of 
effective  properties  in  certain  special  circum¬ 
stances,  the  variational  methods  are  based  on 
and  thus  limited  to  systems  exhibiting  nonlinear 
elastic  or  viscous  material  behavior.  Some  ex¬ 
tensions  have  been  made  to  power-law  materi¬ 
als  and  incremental  theory  of  plasticity.  Also, 
general  agreement  of  the  variational  model  pre¬ 
dictions  has  been  found  with  the  bimodal  yield 
surfaces  of  Section  1.14.3.3  (deBotton,  1995). 
However,  much  more  work  is  needed  before 
more  complex  inelastic  material  behavior, 
which  is  implied  by  available  experiments  dis¬ 


cussed,  for  example,  in  Section  1.14.3.4,  can  be 
represented,  especially  under  a  variable 
loading-unloading  regime. 

Despite  the  progress  witnessed  in  microme¬ 
chanics  of  nonlinear  composite  materials,  there 
are  several  unresolved  issues  that  deserve  future 
attention.  So  far,  constitutive  behavior  of  the 
matrix  has  been  measured  from  specimens  pre¬ 
pared  from  bulk  material.  This  could  be  differ¬ 
ent  from  in  situ  behavior  of  the  matrix,  which  is 
constrained  by  the  reinforcement.  Moreover, 
processing  under  high  temperature  and  pres¬ 
sure  may  cause  a  significant  change  in  the 
local  properties.  Damage,  in  the  form  of  inter¬ 
face  decohesion,  or  microcracks  in  the  phases 
are  often  observed  in  composites  after  fabrica¬ 
tion  or  application  of  service  loads,  and  could 
be  coupled  with  plastic  or  viscoplastic  flow  of 
the  phases.  Finally,  although  periodic  models 
provide  reliable  estimates  of  the  overall  inelas¬ 
tic  properties,  it  may  not  be  representative  of 
certain  composite  systems.  In  this  case,  aver¬ 
aging  models  or  variational  techniques  can  be 
applied  but  in  certain  cases  may  deliver  unreli¬ 
able  predictions.  Macromechanical  models  re¬ 
present  an  alternative,  but  do  not  have  the 
versatility  of  micromechanics. 
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1.14.8  APPENDIX 

Here,  we  record  some  results  related  to 
estimates  of  the  overall  moduli  of  fibrous 
composites  using  averaging  models.  Both  the 
fiber  and  matrix  are  assumed  to  be  transversely 
isotropic  with  a'i  as  the  axis  of  rotational  sym¬ 
metry.  Let  El  and  Er  denote,  respectively,  the 
longitudinal  and  transverse  elastic  Young’s 
modulus,  and  vL  and  vT  denote  Poisson’s 
ratios  under  axial  and  transverse  straining.  Si¬ 
milarly,  we  denote  the  longitudinal  and  trans¬ 
verse  elastic  shear  moduli  by  GL  and  Gy.  The 
elastic  stress-strain  relation  e  =  Mo*  for  either 
phase,  where  M  is  the  elastic  compliance  ma¬ 
trix,  is  written  using  the  engineering  moduli  as 
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The  stiffness  L  =  M  1  is  best  written  in 
terms  of  Hill’s  (1964)  moduli, 

k  =  -[1/Gt  -  4/£x  +  4 vl/Eh] ~\t  =  2fcvL 
n  —  E]^  4-  4&v£  ~  EL-b  f/k,  m  —  GT,/?  =  Cl 

(126) 


The  constitutive  relation  or  =  Le  is  then  writ¬ 
ten  as 


'  crn  ’ 

~n  £  £000" 

£11  ^ 

0*22 

(k  +  m)  (k  —  m)  0  0  0 

£22 

a33 

(k  4  m)  0  0  0 

£33 

a23 

> 

0 

0 

S 

S 

2823 

<*31 

SYM.  p  0 

2e3i 

,  <*12  , 

P. 

.  2£12  . 

(127) 


We  recall  that  Hill’s  formulation  of  the 
Eshelby’s  inclusion  problem,  Section 
1.14.2.2.1,  provides  the  strain  and  stress  con¬ 
centration  factors  in  the  fiber  and  matrix  in 
terms  of  matrices  P  and  Q.  These  matrices 
depend  on  the  shape  of  the  ellipsoidal  inclusion 
and  elastic  moduli  of  the  matrix,  and  are  related 
by  Equation  (20).  For  continuous  fibers  of  a 
circular  cross  section  embedded  in  a  transver¬ 
sely  isotropic  matrix,  the  nonzero  components 
of  P  are  given  by 


coupled  equations  (Hill,  1965a,  1965b), 


Cfkf  ^  cmkm  _ ^  Cffflm  ^mYYlf 

kf  +  m  km+rn  mm  -  m  ny-  -  m 


(129) 


1  _  Cf  Cm  1  _  Cf  Cm 

2 p  p-pm  p-p/k  +  m  kf+m  km  4-  m 

(130) 

Here,  c/  and  cm  are  volume  fractions  of  the 
fiber  and  matrix,  and  mn  pn  and  kr,  r  —  f,m,  are 
Hill’s  moduli  of  the  phases.  The  remaining  two 
moduli,  £  and  /?,  are  found  from  the  universal 
connections 


k  kf  k  km £  cjtf  cm£m  kf  km 

£-£f  £-im  n  -  CfYtf  -  cmnm  if-£m 

(131) 

The  Mori  — Tanaka  model,  on  the  other 
hand,  utilizes  the  form  of  matrix  P  given  in 
Equation  (128)  to  compute  partial  concentra¬ 
tion  factors,  which  provide  the  fiber  average 
stress  and  strain  in  terms  of  their  matrix  coun¬ 
terparts,  Equations  (23)  and  (24).  The  overall 
moduli  are  then  found  from  Equations  (6), 
(23)-(26).  Chen  et  al  (1992)  give  the  Mori- 
Tanaka  moduli  for  fibrous  composites  in  the 
following  explicit  form, 


P22  —  P33  = 


k  1  +  4mi 
%m\(k\  +mi) 


Pl 3  Pn  %mx{k\+m\) 

P44  =  P55  =  —  ,  P<se  = 

2pi 


k\  4-  2m\ 
2m\{k\  +  m\) 


(128) 


2 Cfpmpf  4  Cm  {PmPf  +  Pj) 
2Cfpm  4  Cm  (j?f  4  Pni) 


(132) 


mmrnj{km  4  2  mm)  4  kmmm{cfmf+  cmrnm) 
kmWlm  *4“  {km  4  2?Mm)  {cfY}1m  4  CmtWf) 


(133) 


where  ku  mh  and  px  are  Hill’s  elastic  moduli  of 
the  matrix  phase. 

Estimates  of  the  overall  moduli  of  binary 
composite  materials  by  the  self-consistent 
method  are  found  by  replacing  the  matrix  mod¬ 
uli  in  Equation  (128)  by  the  effective  moduli  of 
the  composite  aggregate,  and  substituting  the 
result  in  Equation  (22).  This  provides  the  con¬ 
centration  factors,  which  are  then  substituted  in 
Equation  (6)  to  determine  the  overall  moduli. 
Numerical  values  of  the  overall  Hill’s  moduli  m, 
p  and  k  are  computed  from  the  following  un- 


j  =  cfkf{km  4  mm)  4  cmkm(kf+rnm) 
Cf{km  4  ftim )  -j-  cm  [kf  4  W/m) 

q _ c/£f{km  4  wi>n)  4*  cm£m (kf  wiffl) 

Cf{km  4  mm)  4  cm  (t kf  4-  mm) 


n  =  cmnm  4-  cjtif  4-  (l  -  cf£f  -  cm£m) 


tf-tm 
kf  -  ktr 


(134) 

(135) 

(136) 


For  completeness,  we  record  the  overall 
Hill’s  moduli  estimated  with  the  vanishing 
fiber  diameter  model  (Dvorak  and  Bahei- 
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El-Din,  1982;  Bahei-El-Din,  1990).  In  this 
model,  the  fiber  constrains  the  matrix  deforma¬ 
tion  in  the  axial  direction  only.  Consequently, 
unit  strain  concentration  factors  are  assumed 
under  axial  straining,  and  unit  stress  concentra¬ 
tion  factors  are  assumed  under  transverse  nor¬ 
mal  and  shear  stresses  as  well  as  longitudinal 
shear  stress.  The  result  is 


=  PmPf 
cfPm  +  Cmpf 


mtnmf 

m  = - - - 

cfmm  +  cmrrif 


k  = 


kmkf 

Cjkm  cmkf 


Cflfkm  "1“  Cmtmkf 
Cfkm  -I-  cmkf 


n  =  cmnm  4*  cfnf  - 


c/cm(ef-em)2 

Cfkm  Cmkj 


(137) 

(138) 


(139) 


